Andrew M. Samuels

20120211

Portfolio

Computational Physics Research

Towards a 27 Day Prediction of Ionospheric and Thermospheric Space Weather

Software/Programming

MATLAB

MATLAB Graphical User Interface (GUI): Projectile Motion

Web Automation

Web Automation with the PHP Programming Language and cURL Module

Twtr - automated microblogging script using Twitter's Rest API. [Download it!]
WPBlogr - automated blogging software, uses SOAP/XML-RPC to post to blogs. [Download it!]
Rubylinkr - a little ditty that automatically creates short URLs
Notebookr - logs into a Google Notebook account, and automatically posts notes

Dynamic, Database Driven Websites with PHP + MySQL

SingleBlackDating.com - a 200,000+ page dating site I made with a huge demographics database
FileFin.com - a freeware file submission website that had 11784 software submissions over its lifetime

Augmented Browsing with Ruby Programming Language + Watir Module

FBPoster
GGrouper
VeohBlack

Python

Java

Various Java Projects #1-12 for Computer Science 111 - Java Programming

Fun Nerdy Stuff

Science Fiction

MATLAB Graphical User Interface (GUI): Projectile Motion

MATLAB is yet another application and programming language I've added to my skill set this year. In this post I'm going to share with you a program I wrote in MATLAB that models projectile motion. It has a convenient little Graphical User Interface (GUI) where the user can input the initial velocity (m/s) and angle (degrees), and the program will spit out a x-y plot.

[Download This!]

20120201

Java Project #12: Statistics


/**
 * class Statistics finds the mean and standard deviation of given numbers.
 * 
 * @author (Andrew Samuels) 
 * @version (04-09-2010)
 */
//import scanner
import java.util.Scanner;

public class Statistics
{
    
    public static void main(String[] args)
    {
        //create variables
        int size;
        int i;
        
        //create scanner
        Scanner keyboard = new Scanner(System.in);
        
        //capture size from user
        System.out.print("How many numbers will be entered: ");
        size = keyboard.nextInt();
        
        //create array of "size"
        double[] p = new double[size];
        
        //read numbers from user one by one
        for(i = 0; i < size; i ++)
        {
            System.out.print("Enter number " + (i + 1) +": ");
            p[i] = keyboard.nextDouble();
        }
        
        //find mean
        double sum = 0;
        for(i = 0; i < size; i++)
        {
            sum = sum + p[i];
        }
        
        double mean = sum/(double)size;
        System.out.println("Mean: " + mean);
        
        //find standard deviation
        double temp = 0;
        for (i = 0; i < size; i++)
        {
            temp = temp + (p[i] - mean)*(p[i] - mean);
        }
        temp = Math.sqrt(temp/size);
        
        System.out.println("Standard deviation: "+ temp);
    }
}

Java Project #11: Array Work

This script shows off various array operations.


/**
 * class ArrayWork.
 * 
 * @author (Andrew Samuels) 
 * @version (04-05-2010)
 */
//import scanner
import java.util.Scanner;

public class ArrayWork
{
    
    public static void main(String[] args)
    {
        //create variables
        int length;
        int i;
        int zeros;
        
        //create scanner
        Scanner keyboard = new Scanner(System.in);
        
        //capture length from user
        System.out.print("Enter the length of array: ");
        length = keyboard.nextInt();
        
        int[] p = new int[length]; //creates array of length, "length"
        System.out.println();
        
        //read integers one by one
        for(i = 0; i < length; i ++)
        {
            System.out.print("Enter number " + (i + 1) +": ");
            p[i] = keyboard.nextInt();
        }
        
        //print numbers in order user entered them
        System.out.print("Numbers: ");
        System.out.println();
        for(i = 0; i < length; i++)
        {
            System.out.print(p[i] + " ");
        }
        
        
        //find sum and average of of numbers and print them
        int sum = 0;
        for(i = 0; i < length; i++)
        {
            sum = sum + p[i];
        }
        System.out.println();
        System.out.println("Sum: " + sum);
        System.out.println();
        
        double average = sum/(double)length;
        System.out.println("Average: " + average);
        System.out.println();
        
        //find max and min in array and print
        int maximum = p[0];
        for (i = 0; i < length; i++)
        {
            if (p[i] > maximum)
            {
                maximum = p[i];
            }
        }
        System.out.println("Maximum: " + maximum);
        System.out.println();
        
        int minimum = p[0];
        for (i = 0; i < length; i++)
        {
            if (p[i] < minimum)
            {
                minimum = p[i];
            }
        }
        System.out.println("Minimum: " + minimum);
        System.out.println();
        
        //print all positive numbers
        
        System.out.println("Positives: ");
        for (i = 0; i < length; i++)
        {
            if (p[i] > 0)
            {
                System.out.print(p[i]+" ");
            }
        }
        
        //find and print zeros
        zeros = countZeros(p);
        System.out.println("Zeros: " + zeros);
        
        //reverse array and print
        printReverse(p);
        
        
    }
    
    //count zeros method
    private static int countZeros(int[] array)
    {
        int answer = 0;
        int i;
        for(i = 0; i < array.length; i++)
        {
            if(array[i] == 0)
            {
                answer = answer + 1;
            }
         
        }
        return answer;
    }
    
    //print reverse method
    private static void printReverse(int[] array)
    {
        int i;
        for(i = array.length - 1; i >= 0; i--)
        {
            System.out.print(array[i] + " ");
        }
        
    }
   
    //display method
    public static void display(int[] array)
    {
        for(int i = 0; i < array.length; i++)
        {
            System.out.println(array[i]);
        }
    }
}

Java Project #10: Area Finder


/**
 * class AreaFinder compute area of square, rectangle, circle.
 * 
 * @author (Andrew Samuels) 
 * @version (03-29-2010)
 */
public class AreaFinder
{
    private static final double pi = 3.1415; //value of constant pi
    private static int count = 0; //count of times area computed
    
    //computes area of square
    public static double squareArea(double sideOfSquare)
    {
        double areaOfSquare = sideOfSquare * sideOfSquare;
        count++;
        return areaOfSquare;
        
    }
    
    //computes area of rectangle
    public static double rectangleArea(double width, double height)
    {
        double areaOfRectangle = width * height;
        count++;
        return areaOfRectangle;
        
    }
    
    //computes area of circle
    public static double circleArea(double radius)
    {
        double areaOfCircle = pi * radius * radius;
        count++;
        return areaOfCircle;
        
    }
    
    //returns number of times the area computing methods were called
    public static int getCallCount()
    {
        return count;
    }
}

20120129

Java Project #9: Pay Raise


/**
 * PayRaise computes pay raise and new salary for employees.
 * 
 * @author (Andrew Samuels) 
 * @version (02-19-2010)
 */

//import scanner
import java.util.Scanner;

public class PayRaise
{

    public static void main(String[] args)
    {
        //create variables
        double currentSalary;
        double payRaise;
        double newSalary;
        String word;
        char performanceRating;
        
        //create scanner
        Scanner keyboard = new Scanner(System.in);
        
        //capture current salary from user
        System.out.print("Please enter current salary: ");
        currentSalary = keyboard.nextDouble();
        
        //capture performance rating from user
        System.out.print("Please enter performance rating (A, B, C): ");
        word = keyboard.next();
        performanceRating = word.charAt(0);
        
        //determine pay raise using switch/case
        switch(performanceRating)
        {
            case 'A':
                payRaise = currentSalary * .06;
                newSalary = currentSalary + payRaise;
                System.out.println("Pay raise is: " + payRaise);
                System.out.println("New salary is: " + newSalary);
                break;
            case 'B':
                payRaise = currentSalary * .04;
                newSalary = currentSalary + payRaise;
                System.out.println("Pay raise is: " + payRaise);
                System.out.println("New salary is: " + newSalary);
                break;
            case 'C':
                payRaise = currentSalary * .02;
                newSalary = currentSalary + payRaise;
                System.out.println("Pay raise is: " + payRaise);
                System.out.println("New salary is: " + newSalary);
                break;
            default:
                System.out.println("Invalid performance rating");
        }
    }
}

Java Project #8 : Tax Program


/**
 * TaxProgram compute tax according to marital status and income.
 * 
 * @author (Andrew Samuels) 
 * @version (02-19-2010)
 */

//import scanner
import java.util.Scanner;

public class TaxProgram
{

    public static void main(String[] args)
    {
        //create variables
        int income;
        int taxAmount;
        String maritalStatus;
        
        //create scanner
        Scanner keyboard =  new Scanner(System.in);
        
        //capture income from user
        System.out.print("Please enter income: ");
        income = keyboard.nextInt();
        
        //capture marital status from user
        System.out.print("Please enter marital status (single, married): ");
        maritalStatus = keyboard.next();
        
        //determine tax amount from marital status and income
        if ((maritalStatus.equals("single") || maritalStatus.equals("Single")) && income < 30000)
        {
            taxAmount = income * 1/5;
            System.out.println("Tax amount is: " + taxAmount);
        }
        else if ((maritalStatus.equals("single") || maritalStatus.equals("Single")) && income > 30000)
        {
            taxAmount = income * 1/4;
            System.out.println("Tax amount is: " + taxAmount);
        }
        else if ((maritalStatus.equals("married") || maritalStatus.equals("Married")) && income < 50000)
        {
            taxAmount = income * 1/10;
            System.out.println("Tax amount is: " + taxAmount);
        }
        else if ((maritalStatus.equals("married") || maritalStatus.equals("Married")) && income > 50000)
        {
            taxAmount = income * 3/20;
            System.out.println("Tax amount is: " + taxAmount);
        }
    }
}

Java Project #7 : Student Grades (Upgraded!)


/**
 * StudentGrades computes the grades of students in a class.
 * 
 * @author (Andrew Samuels) 
 * @version (02-18-2010)
 */

//import scanner
import java.util.Scanner;

public class StudentGrades
{
    
    public static void main(String[] args)
    {
        //create variables
        int homeworkScore;
        int midtermScore;
        int finalexamScore;
        double finalcourseScore;
        
        //create scanner
        Scanner keyboard = new Scanner(System.in);
        
        //capture total homework score from user
        System.out.print("Please enter total homework score (0-200): ");
        homeworkScore = keyboard.nextInt();
        
        //capture midterm score from user
        System.out.print("Please enter midterm score (0-100): ");
        midtermScore = keyboard.nextInt();
        
        //capture final exam score from user
        System.out.print("Please enter final exam score (0-100): ");
        finalexamScore = keyboard.nextInt();
        
        //determine final course score
        finalcourseScore = (50*(homeworkScore/200.0) + 20*(midtermScore/100.0) + 30*(finalexamScore/100.0));
        
        //print final course score
        System.out.println("Final course score is: " + finalcourseScore);
        
        //determine letter grade and meaning of the grade
        if (finalcourseScore >= 90 && finalcourseScore <= 100)
        {
            System.out.println("Letter grade is: A");
            System.out.println("You are: Distinctly above average");
        } 
        else if (finalcourseScore >= 80 && finalcourseScore < 90)
        {
            System.out.println("Letter grade is: B");
            System.out.println("You are: Above average");
        }
        else if (finalcourseScore >= 70 && finalcourseScore < 80)
        {
            System.out.println("Letter grade is: C");
            System.out.println("You are: Average");
        }
        else if (finalcourseScore >= 60 && finalcourseScore < 70)
        {
            System.out.println("Letter grade is: D");
            System.out.println("You are: Below average");
        }
        else if (finalcourseScore >= 0 && finalcourseScore < 60)
        {
            System.out.println("Letter grade is: E");
            System.out.println("You are: Failed");
        }
    }
}

Java Project #6 : Points Distance


/**
 * PointsDistance computes the distance between two points in the two dimensional plane.
 * 
 * @author (Andrew Samuels) 
 * @version (02-04-2010)
 */

//import Scanner from library
import java.util.Scanner;

public class PointsDistance
{
    public static void main(String[] args)
    {
        //create distance formula variables
        double x1;
        double y1;
        double x2;
        double y2;
        double distanceBetweenPoints;
        
        //Create Scanner
        Scanner keyboard = new Scanner(System.in);
        
        //capture x1 from user
        System.out.print("Please enter x coordinate of first point: ");
        x1 = keyboard.nextDouble();
        
        //capture y1 from user
        System.out.print("Please enter y coordinate of first point: ");
        y1 = keyboard.nextDouble();
        
        //capture x2 from user
        System.out.print("Please enter x coordinate of second point: ");
        x2 = keyboard.nextDouble();
        
        //capture y2 from user
        System.out.print("Please enter y coordinate of second point: ");
        y2 = keyboard.nextDouble();
        
        
        //calculate distance between points using distance forumula
        distanceBetweenPoints = Math.sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
        
        
        //print distance between the points
        System.out.println("Distance between (" + x1 + ", " + y1 + ") and (" + x2 + ", " + y2 +") is: " + distanceBetweenPoints);
    }
}

Java Project #5 : Bank Investment


/**
 * BankInvestment computes the total money earned in an investment bank.
 * 
 * @author (Andrew Samuels) 
 * @version (02-04-2010)
 */
//importing Scanner from library
import java.util.Scanner;

public class BankInvestment
{
    public static void main(String[] args)
    {
    //create variables
    int initialInvestmentMoney; //initial investment money in dollars
    double interestRate; //interest rate in percentage
    int numberOfYears; //number of years of investment
    double totalMoneyEarnedDouble; //total money earned, double type, (converted to integer later)
    
    //Create Scanner
    Scanner keyboard = new Scanner(System.in);
        
    //capture initial investment from user
    System.out.print("Please enter initial investment (dollars): ");
    initialInvestmentMoney = keyboard.nextInt();
    
    //capture interest rate from user
    System.out.print("Please enter interest rate (percent): ");
    interestRate = keyboard.nextDouble();
    
    //capture numbers of years of investment from user
    System.out.print("Please enter years of investment (years): ");
    numberOfYears = keyboard.nextInt();
    
    //Calculate total money earned with formula p(1 + r/100)^n
    totalMoneyEarnedDouble = initialInvestmentMoney*(Math.pow(1 + interestRate/100, numberOfYears));
    
    //convert total money earned from double to integer
    int totalMoneyEarned = (int)totalMoneyEarnedDouble;
    
    //print total money earned
    System.out.println("The total money earned (dollars): " + totalMoneyEarned);
    }
}

20120128

Java Project #4: Travel Cost

/**
 * TravelCost computes the cost of travel using a vehicle.
 * 
 * @author (Andrew Samuels) 
 * @version (02-04-2010)
 */

//importing Scanner from library
import java.util.Scanner;

public class TravelCost
{
    public static void main(String[] args)
    {
    //create variables
    int travelDistance; //travel distance in miles
    int gasolineUsage; //gasoline usage in miles per gallon
    double gasolinePrice; //gasoline price in dollars per gallon
    double totalTravelCost; //The total cost of traveling
    
    //Create Scanner
    Scanner keyboard = new Scanner(System.in);
        
    //capture travel distance from user
    System.out.print("Please enter travel distance (miles): ");
    travelDistance = keyboard.nextInt();
    
    //capture gas usage from user
    System.out.print("Please enter gas usage of vehicle (miles per gallon): ");
    gasolineUsage = keyboard.nextInt();
    
    //capture gas price from user
    System.out.print("Please enter gas price (dollars per gallon): ");
    gasolinePrice = keyboard.nextDouble();
    
    //Calculate
    totalTravelCost = (travelDistance/gasolineUsage)*gasolinePrice;
    
    //print total travel cost
    System.out.println("The total travel cost (dollars): " + totalTravelCost);
    }
}

Java Project #3 : Student Grade

/**
 * StudentGrade determines the final semester score of students in a course.
 * 
 * @author (Andrew Samuels) 
 * @version (02-01-2010)
 */
//import scanner class
import java.util.Scanner;

public class StudentGrade
{

    public static void main(String[] args)
    {
        //initialize variables
        int programs;
        int homeworks;
        int examinations;
        double finalSemesterScore;
        
        //create scanner
        Scanner keyboard = new Scanner(System.in);
        
        //input total program score from the user
        System.out.print("Please enter total program score (0-200): ");
        programs = keyboard.nextInt();
        
        //input total homework score from the user
        System.out.print("Please enter total homework score (0-500): ");
        homeworks = keyboard.nextInt();
        
        //input examination score from the user
        System.out.print("Please enter examination score (0-100): ");
        examinations = keyboard.nextInt();
        
        //calculate final semester score
        finalSemesterScore = 50*(programs*.005) + 30*(homeworks*.002) + 20*(examinations*.01);
        
        //print final semester score
        System.out.println("Final semester score is: " + finalSemesterScore);
        
        
    }
}

Java Project #2 : Coinjar

Got spare change in the piggie jar? Coinjar is a Java script that determines the value of coins in a jar. Tell it how many quarters, dimes, nickels, and pennies you have and it will spit out the dollar value.

/**
 * CoinJar determines the value of coins in a jar.
 * 
 * @author (Andrew Samuels) 
 * @version (02-01-2010)
 */
//import scanner class
import java.util.Scanner;

public class CoinJar
{

    public static void main(String[] args)
    {
        //create currency variables in pennies
        int quarters = 25;
        int dimes = 10 ;
        int nickels = 5;
        int pennies = 1;
        
        //create variables for user number of change
        int quartersNumber;
        int dimesNumber;
        int nickelsNumber;
        int penniesNumber;
        
        //create variables for value of change
        int quartersTotal;
        int dimesTotal;
        int nickelsTotal;
        int penniesTotal;
        
        //variables for total value as double
        int totalValue;
        double dollarCents;
        
        
        //create scanner
        Scanner keyboard = new Scanner(System.in);
        
        //input quarters from the user
        System.out.print("Please enter number of quarters: ");
        quartersNumber = keyboard.nextInt();
        
        //input dimes from the user
        System.out.print("Please enter number of dimes: ");
        dimesNumber = keyboard.nextInt();
        
        //input nickels from the user
        System.out.print("Please enter number of nickels: ");
        nickelsNumber = keyboard.nextInt();
        
         //input pennies from the user
        System.out.print("Please enter number of pennies: ");
        penniesNumber = keyboard.nextInt();
        
        //calculate quarters value
        quartersTotal = quarters*quartersNumber;
    
        //calculate dimes value
        dimesTotal = dimes*dimesNumber;
        
        //calculate nickels value
        nickelsTotal = nickels*nickelsNumber;
        
        //calculate pennies value
        penniesTotal = pennies*penniesNumber;
        
        //calculate total value in pennies
        totalValue = quartersTotal + dimesTotal + nickelsTotal + penniesTotal;
        
        //calculate total value as double in dollarCents ex. 2.19
        dollarCents = totalValue*(.01);
        
        //convert dollarCents to String for later use.
        String dollarCentsStr;
        dollarCentsStr = "" + dollarCents;
        
        //convert to dollarCents into integer
        int dollarCentsInt = (int)dollarCents;
        
        //separate cents using substring
        String number = dollarCentsStr;
        String result;
        result = number.substring(2);
       
        //print total dollars and cents
        System.out.println("Jar has " + dollarCentsInt + " dollars and " + result + " cents");
        
        
    }
}

Java Project #1 : The Sphere Problem

This is a straightforward problem. This little ditty computes the volume and surface area of a sphere. Java is fun, no?

/**
 * SphereProblem reads the radius of a sphere from the user and computes its volume and surface area, then prints them out..
 * 
 * @author (Andrew Samuels) 
 * @version (02-01-2010)
 */
//import scanner class
import java.util.Scanner;

public class SphereProblem
{

    public static void main(String[] args)
    {
        //initialize variables
        double radius; //radius of a sphere
        double volume; //volume of a sphere
        double area; //area of a sphere
        double PI = 3.1415;
        
        //create scanner
        Scanner keyboard = new Scanner(System.in);
        
        //input radius from the user
        System.out.print("Please enter radius of a sphere: ");
        radius = keyboard.nextDouble();
        
        //calculate volume
        volume = (4/3)*PI*radius*radius*radius;
    
        //calculate area
        area = 4*PI*radius*radius;
        
        //print volume
        System.out.println("Volume of sphere is: " + volume);
        
        //print area
        System.out.println("Area of sphere is: " + area);
        
    }
}

20111228

Fun with Python: Simulating Baseball Pitches using 4th Order Runge-Kutta Method

Here is a description of the problem in PDF format. Basically this code simulates four different baseball pitches (slider, curve, fastball, screwball) by solving differential equations using the 4th order Runge-Kutta method.

#!/opt/local/bin/python
#Andrew Samuels
#baseball.py

from scipy import *
from pdb import *


#create arrays
dt = 1e-4
nMAX = int(1/dt)

x = zeros(nMAX)
y = zeros(nMAX)
z = zeros(nMAX)
vx = zeros(nMAX)
vy = zeros(nMAX)
vz = zeros(nMAX)
v = zeros(nMAX)


#constants
B = 4.1e-4 #magnus force (dimensionless)
w = 1800/60.0 * 2*pi #rotaton rate of baseball, 1800 rpm to rad/s
g = 9.81 #acceleration of gravity (m/s^2)
v0_fb = 42.4688 #initial speed of pitch, fastball (m/s)
v0_op = 38.0 #initial speed of pitch, other pitches (m/s)
l = 18.44 #distance from pitcher (m)
T_fb = l/v0_fb #step-length, estimated flight time, fastball
T_op = l/v0_op  #step-length, estimated flight time, other pitches
nT_fb = int(T_fb/dt) #number of iterations for fastball
nT_op = int(T_op/dt) #number of iterations for other pitches
d2r = pi/180 # convert degrees to radians
theta = 1 * d2r # elevation angle of pitch, degrees
phi1 = 225 * d2r #direction of w, relative to z, fastball
phi2 = 45 * d2r #direction of w, relative to z, curveball
phi3 = 0 * d2r #direction of w, relative to z, slider
phi4 = 135 * d2r #direction of w, relative to z, screwball
#drag force constants
c1 = 0.0039
c2 = 0.0058
vd = 35.0 #(m/s)
delta = 5.0 # (m/s)

#initial conditions
x[0] = 0.0
y[0] = 0.0
z[0] = 0.0
vx[0] = v0_fb * cos(theta) 
vy[0] = 0.0
vz[0] = v0_fb * sin(theta)


#pass dx/dt=vx through runkut4
def rk4dx(vx,dt):
    def dxdt(vx):
        dxdt = vx
        return dxdt
    K0 = dt * dxdt(vx)
    K1 = dt * dxdt(vx + K0/2.0)
    K2 = dt * dxdt(vx + K1/2.0)
    K3 = dt * dxdt(vx + K2)
    vxnew = (K0 + 2.0*K1 + 2.0*K2 + K3)/6.0
    return vxnew

#pass dy/dt=vy through runkut4
def rk4dy(vy,dt):
    def dydt(vy):
        dydt = vy
        return dydt
    K0 = dt * dydt(vy)
    K1 = dt * dydt(vy + K0/2.0)
    K2 = dt * dydt(vy + K1/2.0)
    K3 = dt * dydt(vy + K2)
    vynew = (K0 + 2.0*K1 + 2.0*K2 + K3)/6.0
    return vynew
    
#pass dz/dt=vz through runkut4
def rk4dz(vz,dt):
    def dzdt(vz):
        dzdt = vz
        return dzdt
    K0 = dt * dzdt(vz)
    K1 = dt * dzdt(vz + K0/2.0)
    K2 = dt * dzdt(vz + K1/2.0)
    K3 = dt * dzdt(vz + K2)
    vznew = (K0 + 2.0*K1 + 2.0*K2 + K3)/6.0
    return vznew

#solve accelerations with runge kutta 4
def rk4(v,vx,vy,vz,dt,phi):
    def ax(v,vx,vy,vz,phi):
    #drag force function, Fd(v)
        def Fd(v):
            return c1 + (c2/(1.0 + exp((v-vd)/delta)))
        return -Fd(v)*v*vx + B*w*(vz*sin(phi)-vy*cos(phi))
    def ay(v,vx,vy,phi):
    #drag force function, Fd(v)
        def Fd(v):
            return c1 + (c2/(1.0 + exp((v-vd)/delta)))
        return -Fd(v)*v*vy + B*w*vx*cos(phi1)
    def az(v,vx,vz,phi):
    #drag force function, Fd(v)
        def Fd(v):
            return c1 + (c2/(1.0 + exp((v-vd)/delta)))
        return -g*-Fd(v)*v*vz - B*w*vx*sin(phi)
    K0 = dt * ax(v,vx,vy,vz,phi)
    L0 = dt * ay(v,vx,vy,phi)
    M0 = dt * az(v,vx,vz,phi)
    K1 = dt * ax(v + dt/2.0,vx + K0/2.0,vy + L0/2.0,vz + M0/2.0,phi)
    L1 = dt * ay(v + dt/2.0,vx + K0/2.0,vy + L0/2.0,phi)
    M1 = dt * az(v + dt/2.0,vx + K0/2.0,vz + M0/2.0,phi)
    K2 = dt * ax(v + dt/2.0,vx + K1/2.0,vy + L1/2.0,vz + M1/2.0,phi)
    L2 = dt * ay(v + dt/2.0,vx + K1/2.0,vy + L1/2.0,phi)
    M2 = dt * az(v + dt/2.0,vx + K1/2.0,vz + M1/2.0,phi)
    K3 = dt * ax(v + dt,vx + K2,vy + L2,vz + M0,phi)
    L3 = dt * ay(v + dt,vx + K1,vy + L2,phi)
    M3 = dt * az(v + dt,vx + K1,vz + M2,phi)
    axnew = (K0 + 2.0*K1 + 2.0*K2 + K3)/6.0
    aynew = (L0 + 2.0*L1 + 2.0*L2 + L3)/6.0
    aznew = (M0 + 2.0*M1 + 2.0*M2 + M3)/6.0
    return axnew,aynew,aznew


f=open('fastball.dat','w')

#first time step is n + 1
#start loop

for i in range(0,nT_fb):
    v = (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i])
    v = (v)**0.5
    print "X: " + str(x[i]) + " , Y: " + str(y[i]) + " , Z: " + str(z[i]) + " ,T: " + str(x[i]/v) +" , V: "+ str(v) + "\n"
    f.write(str(x[i])+'\t'+str(y[i])+'\t'+ str(z[i])+'\t'+str(x[i]/v) + '\t' + str(v) +'\n')
    x[i+1] = x[i] + rk4dx(vx[i],dt) 
    y[i+1] = y[i] + rk4dy(vy[i],dt)
    z[i+1] = z[i] + rk4dz(vz[i],dt)
    vx[i+1] = vx[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi1)[0]
    vy[i+1] = vy[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi1)[1]
    vz[i+1] = vz[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi1)[2]


f.close()
print "Created fastball.dat"

#change initial conditions for the other pitches
vx[0] = v0_op * cos(theta) 
vz[0] = v0_op * sin(theta)

#curveball
h=open('curveball.dat','w')

for i in range(0,nT_op):
    v = (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i])
    v = (v)**0.5
    print "X: " + str(x[i]) + " , Y: " + str(y[i]) + " , Z: " + str(z[i]) + " ,T: " + str(x[i]/v) +" , V: "+ str(v) + "\n"
    h.write(str(x[i])+'\t'+str(y[i])+'\t'+ str(z[i])+'\t'+str(x[i]/v) + '\t' + str(v) +'\n')
    x[i+1] = x[i] + rk4dx(vx[i],dt) 
    y[i+1] = y[i] + rk4dy(vy[i],dt)
    z[i+1] = z[i] + rk4dz(vz[i],dt)
    vx[i+1] = vx[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi2)[0]
    vy[i+1] = vy[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi2)[1]
    vz[i+1] = vz[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi2)[2]


h.close()
print "Created curveball.dat"

#slider
j=open('slider.dat','w')

for i in range(0,nT_op):
    v = (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i])
    v = (v)**0.5
    print "X: " + str(x[i]) + " , Y: " + str(y[i]) + " , Z: " + str(z[i]) + " ,T: " + str(x[i]/v) +" , V: "+ str(v) + "\n"
    j.write(str(x[i])+'\t'+str(y[i])+'\t'+ str(z[i])+'\t'+str(x[i]/v) + '\t' + str(v) +'\n')
    x[i+1] = x[i] + rk4dx(vx[i],dt) 
    y[i+1] = y[i] + rk4dy(vy[i],dt)
    z[i+1] = z[i] + rk4dz(vz[i],dt)
    vx[i+1] = vx[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi3)[0]
    vy[i+1] = vy[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi3)[1]
    vz[i+1] = vz[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi3)[2]


j.close()
print "Created slider.dat"

#screwball
o=open('screwball.dat','w')

for i in range(0,nT_op):
    v = (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i])
    v = (v)**0.5
    print "X: " + str(x[i]) + " , Y: " + str(y[i]) + " , Z: " + str(z[i]) + " ,T: " + str(x[i]/v) +" , V: "+ str(v) + "\n"
    o.write(str(x[i])+'\t'+str(y[i])+'\t'+ str(z[i])+'\t'+str(x[i]/v) + '\t' + str(v) +'\n')
    x[i+1] = x[i] + rk4dx(vx[i],dt) 
    y[i+1] = y[i] + rk4dy(vy[i],dt)
    z[i+1] = z[i] + rk4dz(vz[i],dt)
    vx[i+1] = vx[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi4)[0]
    vy[i+1] = vy[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi4)[1]
    vz[i+1] = vz[i] + rk4(v,vx[i],vy[i],vz[i],dt,phi4)[2]


o.close()
print "Created screwball.dat"

Fun with Python: Electric Field Using Trapezoid Method

Suppose a line of charge is given that lies on the y-axis from y=-10 cm to y=+10 cm. At any point, the electric field can be calculated by the following function:

E_x=(sigma/(4*pi*epsilon_0))*intregral_-10_10(x/((x^2+y^2)^(3/2)))dy

where sigma is the linear charge density.

This python program solves this equation using the trapezoid method for numerical integration and outputs the data into a .dat file.

#Andrew Samuels
#efield.py
#find the electric field for a given function
#Uses the trapezoid method

#import modules
from math import *

#open a data file
dfile = open('efdata.dat','w')

#make some functions
def efunc(x,y):
    return  x/(x**2.+y**2.)**(3/2.)

def efTrap(f,a,b,step=1):
    A = -(f(x,a) + f(x,b))
    S = float(b - a)
    for n in range(step+1):
        y = a + n*S/step
        A += 2*f(x,y)
    return A*float(b-a)/(2*step)

#constant
k = 10.5/1.113e-4 #10.5microC/4*pi*e_0*1e6

#integrate from x=1 to 100 in steps of 5
x=1
while x <= 100 :
    ef = k * efTrap(efunc,-10,10,10000)
    print "x = " + str(x) + " , eField = " + str(ef)
    dfile.write(str(x) + '\t' + str(ef) + '\n')
    x += 5

#close data file
dfile.close()
print "Done! Created efdata.dat"

Fun with Python: Linear Algebra Solver

This program is similar to my Gaussian Elimination script only this time the matrix is solved using Back Substitution, then checked for accuracy using the Sci-Py module.

#Andrew Samuels
#linear equation solver. Does Gauss-Jordon elimination, then back substituion.

from scipy import *
from scipy import linalg
from pdb import *

#function to compute dot product
def dot(a, b):
    return reduce(lambda sum, p: sum + p[0]*p[1], zip(a,b), 0)

matrix = zeros((3,4))
#aug = zeros((3,1))

matrix[0,0] = 1
matrix[0,1] = 1
matrix[0,2] = 1
matrix[0,3] = 0
matrix[1,0] = 1
matrix[1,1] = -2
matrix[1,2] = 2
matrix[1,3] = 4
matrix[2,0] = 1
matrix[2,1] = 2
matrix[2,2] = -1
matrix[2,3] = -1


pivot = 0
ncols = 4
nrows = 3


#for all pivots
for col in range(ncols-1):
    #divide pivot by itself to make it 1
    matrix[pivot,:] = matrix[pivot,:]/matrix[pivot,pivot]
   
    #for all pivots that haven't been done yet
    for row in range(pivot+1,nrows):
        matrix[row,:] = matrix[row,:] - matrix[pivot,:]
    pivot = pivot + 1

print " "
print "Row Echelon Form: "   
print matrix 
print "\n"

#solve the matrix and print solution
for pivot in range(nrows-1,-1,-1):
    matrix[pivot] = (matrix[pivot] - dot(matrix[pivot,pivot+1:nrows],matrix[pivot+1:nrows]))/matrix[pivot,pivot]
print "Solved by Back Subst.: "
print matrix 
print "\n"  

#check with scipy
a = array([[1,1,1],[1,-2,2],[1,2,-1]])
b = array([[0],[4],[-1]])
sol = linalg.solve(a,b)
print "Solved by Scipy (check): "
print sol
print "\n"

Fun with Python: Gaussian Elimination

Here is my python program that does Gaussian Elimination on a 3 x 4 matrix. Of course, there are existing modules with functions that could do this automatically, but this was a good learning exercise doing it the hard way.

#Andrew Samuels
#Gaussian Elimination
#performs gaussian elimination on a 3 by 4 matrix

#import scipy module
from scipy import *
from pdb import *

#construct the matrix
matrix = zeros((3,4))

matrix[0,0] = 1
matrix[0,1] = 1
matrix[0,2] = 1
matrix[0,3] = 0
matrix[1,0] = 1
matrix[1,1] = -2
matrix[1,2] = 2
matrix[1,3] = 4
matrix[2,0] = 1
matrix[2,1] = 2
matrix[2,2] = -1
matrix[2,3] = -1


pivot = 0
ncols = 4
nrows = 3

#for all pivots
for col in range(ncols-1):
    #divide pivot by itself to make it 1
    matrix[pivot,:] = matrix[pivot,:]/matrix[pivot,pivot]
   
    #for all pivots that haven't been done yet
    for row in range(pivot+1,nrows):
        matrix[row,:] = matrix[row,:] - matrix[pivot,:]
    pivot = pivot + 1
   

print matrix

Fun with Python: Factorial Calculator

This program written in python calculates a factorial when given an integer. The program will complain if you don't give it a valid integer, such as a letter or dollar sign.

#Andrew Samuels
#factorial calculator
#Fall 2011
#calculates a factorial from integer input from user

n = raw_input('Enter an integer: ')
try:
    n=int(n)
except:
    print 'Enter a valid integer, dummy!'
    exit(1)

f=1
while n > 1 :
    f = f * n
    n = n - 1
    print f

20110819

Towards a 27 Day Prediction of Ionospheric and Thermospheric Space Weather

In March of 2011, during my winter semester at Eastern Michigan University Department of Physics & Astronomy, I had the honor of writing a proposal for a research scholarship program award with my computational physics professor, Dr. Dave Pawlowski. The proposal was for the University Research Stimulus Program (URSP) at EMU. The objective of this award was to facilitate research partnerships between undergraduate students and Eastern Michigan University faculty along with a reward of scholarship for the undergraduate.

Shortly after submitting my proposal, I had received confirmation that I had in fact been granted the URSP award. During the months of April through August 2011, I had been working daily as an undergraduate research assistant under the direction of Dr. Pawlowski. This research was my first plunge into the discipline of computational physics. Over this time period I had gained experience in the subject and performed the following:

• Running computer simulations of Earth’s upper atmosphere based on observations of the sun. Creating a probabilistic forecast for space weather.
• Programmed software scripts using IDL and FORTRAN
• Completed NASA’s Information Security Training, granted access to NASA’s Pleiades supercomputer.
• Parallel computing protocol with Message Passing Interface (MPI)

This post is the contents of a research poster I had put together for a presentation I had made to three groups of visiting 7th and 8th graders, fellow science undergraduates and faculty. You may also download my poster in PDF format here.

Abstract

The prospect of what we can do to increase the knowledge of space weather is promising. One of the practical reasons for doing so is the fact that perturbations in density and temperature in the upper atmosphere can have a significant effect on satellites and instrumentation in low orbit of Earth. By utilizing ensembles of simulations or ensembles that are based on the uncertainty of the system drivers, we can run computer simulations to help us make probabilistic forecasts about space weather. In this study, we prepared an ensemble simulation based on uncertainties in the solar extreme ultraviolet (EUV) flux. Then, using the Global Ionosphere-Thermosphere Model (GITM) we ran simulations simultaneously using these ensembles, and analyzed the results.

Introduction

•Space weather

•the conditions in Earth’s upper atmosphere (ionosphere and thermosphere) in response to changing conditions on the sun.

Why do we scarcely hear meteorologists in our local news reports making predictions about space weather?

•The overarching reasons:

•Society at large generally doesn’t know that space weather affects expensive infrastructure we use daily (e.g. satellite communications, GPS navigation)

•Space weather forecasting is currently in significant need of more scientific research in order to make reasonably accurate, long-term predictions.

•Ensemble forecasting

•Create a statistical sample of weather outcomes by performing multiple computer simulations simultaneously.

•Ensembles, are based on a number of input conditions called “system drivers”.

•Ensemble members are created based on the uncertainty of these drivers and are used as inputs to our Global Ionosphere and Thermosphere Model (GITM) and simulations are performed which will create a probabilistic forecast of space weather.

This first segment of this long-term study involves understanding and collecting solar flux data. Utilizing the known uncertainties of this data, we created an ensemble of drivers that were used as input to multiple GITM simulations. This study examines the effect of the uncertainty in the solar flux on the state of the upper atmosphere.

About the Global Ionosphere and Thermosphere Models

•3-D coupled thermosphere-ionosphere model [Ridley et al., 2006]

•Solves for 11 Neutral and 10 Ion species, neutral winds (horizontal and vertical), ion and electron velocities and neutral, ion, and electron temperatures

•Does not assume hydrostatic equilibrium: Coriolis, vertical ion drag, non- constant gravity, massive auroral zone heating

•Flexible grid resolution, fully parallel, with an altitude grid

•Allows for 1D simulations to perform long duration

For this study: 2.5o Latitude x 5o Longitude resolution using 64 processors on the NASA Pleiades supercomputer

Drivers: Solar Flux

Plot of the Solar Flux vs. Time. The black line represents the actual solar flux. The colored lines represent the 5 ensemble members that are used to drive the simulations.

Thermosphere Results

The following plots show the mean mass atmospheric density which is determined from the mean of 5 ensemble members (top plots) and their standard deviations (bottom plots) of the thermosphere at a 403 km altitude.

1:00 UT: The simulation has just started and thus there is small difference between the simulations due to the different solar fluxes being used.

4:00 UT: The effects of the flare on the mass density begin to appear.

6:00 UT: Here the effects of the flare reaches its maximum. On a large scale, this is when the largest effects occur. Notice the time delay in effects from Figure 1's peak.

13:00 UT: The effects of the uncertainty have propagated globally. High levels of standard deviation are probably due to significantly different wind patterns. The black arrow points to a maximum uncertainty of 1.7%

Ionosphere Results

The following plots show the mean electron density which is determined from the mean of the 5 ensemble members (top plots) and their standard deviations (bottom plots) of the ionosphere at an altitude of 118 km (left) and 403km (right) .

3:00 UT (118 km): Flare has just ended. Ionosphere is significantly uncertain on the day side.

3:00 UT (403 km): At this higher altitude our standard deviation has significantly increased along with the mean density of electrons.

Conclusions and Future Work

Given the uncertainties in the solar flux, it takes several hours for thermosphere to become uncertain. Within 12 hours, the uncertainties have propagated globally in the thermosphere. This is interesting because the thermosphere only affects the day side directly.

•At 13 UT the maximum uncertainty is 1.7%

The ionosphere, however can become uncertain very quickly.

•At 3:00 UT (403 km): The maximum uncertainty is 47%.

The ultimate goal of this research is to predict space weather over one solar rotation period (27 days). We need to study the uncertainty involved in using data on this month’s drivers to predict next month’s space weather. This research also has yet to account for both types of uncertainty in the driver: flare and daily uncertainty. In the future, it will use even more ensemble inputs, like interplanetary magnetic fields (IMF) and solar wind to the GITM model.

References and Acknowledgements

Ridley, A.J., Deng, Y., Toth, G., The Global Ionosphere-Thermosphere Model, J. Atmos. Sol-Terr. Phys., 68, 839, 2006

Flare irradiance data was processed and provided by Phil Chamberlin and Anne Wilson, University of Colorado. http://lasp.colorado.edu/lisird/fism/fism.html

Special thanks to Eastern Michigan University for proving the Undergraduate Research Stimulus Program.