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. 
/**
* 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);
}
}
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]);
}
}
}
/**
* 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;
}
}
/**
* 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");
}
}
}
/**
* 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);
}
}
}
/**
* 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");
}
}
}
/**
* 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);
}
}
/**
* 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);
}
}
/**
* 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);
}
}
/**
* 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);
}
}
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");
}
}
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);
}
}
#!/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"
#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"
#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"
#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
#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