#!/usr/bin/env python3
# -*- coding: utf-8 -*-

###
# Name: Nikki and Riley
# Student ID: 2267883 2274503
# Email: schwa218@mail.chapman.edu keda106@mail.chapman.edu
# Course: PHYS220/MATH220/CPSC220 Fall 2017
# Assignment: Classwork 4
###

"""Classwork 04
This classwork introduces numpy arrays and compares their performance to
python lists.
"""

import math
import numpy as np

#GAUSSIAN LIST

def gen_gaussian_list(a, b, n=1000):
    """gen_gaussian_list(a, b, n=1000)
    Generate a discrete approximation of a Gaussian function, including its
    domain and range, stored as a pair of vanilla python lists.
    
    Args:
        a (float) : Lower bound of domain
        b (float) : Upper bound of domain
        n (int, optional) : Number of points in domain, defaults to 1000.
    
    Returns:
        (x, g) : Pair of lists of floats
            x  : [a, ..., b] List of n equally spaced floats between a and b
            g  : [g(a), ..., g(b)] List of Gaussian values matched to x
    """
    dx = (b-a)/(n-1)                         # spacing between points
    x = [a + k*dx for k in range(n)]         # domain list
    
    # Local implementation of a Gaussian function
    def gauss(x):
        return (1/math.sqrt(2*math.pi))*math.exp(-x**2/2)
    
    g = [gauss(xk) for xk in x]                  # range list
    return (x, g)

#NEW GAUSS ARRAY USING NUMPY

def gen_gaussian_array(a, b, n=1000):
    """gen_gaussian_array(a, b, n=1000)
    Generate a discrete approximation of a Gaussian function, including its
    domain and range, stored as a pair of numpy arrays.
    
    Args:
        a (float) : Lower bound of domain
        b (float) : Upper bound of domain
        n (int, optional) : Number of points in domain, defaults to 1000.
    
    Returns:
        (x, g) : Pair of numpy arrays of float64
            x  : [a, ..., b] Array of n equally spaced float64 between a and b
            g  : [g(a), ..., g(b)] Array of Gaussian values matched to x
    """
    ### No need for np.array here. np.linspace already returns an array.
    x = np.array(np.linspace(a,b,n),dtype=np.float64)
    def gauss(x):
        return (1/math.sqrt(2*math.pi))*math.exp(-x**2/2)  ### Do not use math if using numpy. Use np.sqrt etc.
    ### Vectorize is not needed here since all operations are already vectorized (if np.exp and np.pi and np.sqrt are used)
    fx = np.vectorize(gauss)
    g = fx(x)  ### No need to redefine g here. Just return gauss(x) directly.
    return(x,g)

#SINC LIST

def sinc_list(a, b, n=10000):
    """sinc_list(a, b, n=10000)
    Generate a discrete approximation of a Sinc function, including its
    domain and range, stored as a pair of vanilla python lists.
    
    Args:
        a (float) : Lower bound of domain
        b (float) : Upper bound of domain
        n (int, optional) : Number of points in domain, defaults to 10000.
    
    Returns:
        (x, g) : Pair of lists of floats
            x  : [a, ..., b] List of n equally spaced floats between a and b
            g  : [g(a), ..., g(b)] List of Gaussian values matched to x
    """
    dx = (b-a)/(n-1)                         # spacing between points
    x = [a + k*dx for k in range(n)]         # domain list
    
    # Local implementation of a Sinc function
    def sinc(x):
        if x == 0:
           return 1
        else:
            return (math.sin(x)/x)
    
    g = [sinc(xk) for xk in x]                  # range list
    return (x, g)


#NEW SINC ARRAY USING NUMPY

def sinc_array(a, b, n=10000):
    """sinc_array(a, b, n=10000)
    Generate a discrete approximation of a Sinc function, including its
    domain and range, stored as a pair of numpy arrays.
    
    Args:
        a (float) : Lower bound of domain
        b (float) : Upper bound of domain
        n (int, optional) : Number of points in domain, defaults to 10000.

    Returns:
        (x, g) : Pair of numpy arrays of float64
            x  : [a, ..., b] Array of n equally spaced float64 between a and b
            g  : [g(a), ..., g(b)] Array of Gaussian values matched to x
    """
    x = np.array(np.linspace(a,b,n),dtype=np.float64)  ### no need for np.array
    def sinc(x):
        if x == 0:
            return 1
        else:
            return (math.sin(x)/x)
    fx = np.vectorize(sinc)  ### Vectorize IS needed here, since the if statement in sinc is not vectorized yet.
    g = fx(x)  ### no need for g, just return fx(x) directly.
    return(x,g)


def main(a,b,n=1000):
    """main(a, b, n=1000)
    Main function for command line operation. Prints result of Gaussian to screen.
    
    Args:
        a (float) : Lower bound of domain
        b (float) : Upper bound of domain
        n (int, optional) : Number of points in domain, defaults to 1000.
    
    Returns:
        None
    
    Effects:
        Prints Gaussian to screen.
    """
    # You can unpack tuples as return values easily
    x, g = gen_gaussian_list(a,b,n)
    # The zip function takes two lists and generates a list of matched pairs
    for (xk, gk) in zip(x, g):
        # The format command replaces each {} with the value of a variable
        print("({}, {})".format(xk, gk))


# Protected main block for command line operation
if __name__ == "__main__":
    # The sys module contains features for running programs
    import sys
    # The sys.argv list variable contains all command line arguments
    #    sys.argv[0] is the program name always
    #    sys.argv[1] is the first command line argument, etc
    if len(sys.argv) == 4:
        a = float(sys.argv[1])
        b = float(sys.argv[2])
        n = int(sys.argv[3])
        main(a,b,n)
    elif len(sys.argv) == 3:
        a = float(sys.argv[1])
        b = float(sys.argv[2])
        main(a,b)
    else:
        print("Usage: cw04.py a b [n]")
        print("  a : float, lower bound of domain")
        print("  b : float, upper bound of domain")
        print("  n : integer, number of points in domain")
        exit(1)

