Personel Information

BSc: Physics Department, Middle East Technical University (METU), Ankara, Turkey
MSc: Physics Department, METU
PhD: Qualified, Physics Department, METU


Monday, July 11, 2016

Second Harmonic Generation N=1:21



gnuplot> set xrange [-180:180]
gnuplot> set yrange [-180:180]


splot sin(cos(x*pi/180))*sin(cos(x*pi/180))/(cos(x*pi/180)*cos(x*pi/180))*sin(cos(y*pi/180))*sin(cos(y*pi/180))/(cos(y*pi/180)*cos(y*pi/180))+sin(cos((x+1)*pi/180))*sin(cos((x+1)*pi/180))/(cos((x+1)*pi/180)*cos((x+1)*pi/180))*sin(cos((y+1)*pi/180))*sin(cos((y+1)*pi/180))/(cos((y+1)*pi/180)*cos((y+1)*pi/180))+sin(cos((x+2)*pi/180))*sin(cos((x+2)*pi/180))/(cos((x+2)*pi/180)*cos((x+2)*pi/180))*sin(cos((y+2)*pi/180))*sin(cos((y+2)*pi/180))/(cos((y+2)*pi/180)*cos((y+2)*pi/180))+sin(cos((x+3)*pi/180))*sin(cos((x+3)*pi/180))/(cos((x+3)*pi/180)*cos((x+3)*pi/180))*sin(cos((y+3)*pi/180))*sin(cos((y+3)*pi/180))/(cos((y+3)*pi/180)*cos((y+3)*pi/180))+sin(cos((x+4)*pi/180))*sin(cos((x+4)*pi/180))/(cos((x+4)*pi/180)*cos((x+4)*pi/180))*sin(cos((y+4)*pi/180))*sin(cos((y+4)*pi/180))/(cos((y+4)*pi/180)*cos((y+4)*pi/180)) +sin(cos((x+5)*pi/180))*sin(cos((x+5)*pi/180))/(cos((x+5)*pi/180)*cos((x+5)*pi/180))*sin(cos((y+5)*pi/180))*sin(cos((y+5)*pi/180))/(cos((y+5)*pi/180)*cos((y+5)*pi/180))+sin(cos((x+6)*pi/180))*sin(cos((x+6)*pi/180))/(cos((x+6)*pi/180)*cos((x+6)*pi/180))*sin(cos((y+6)*pi/180))*sin(cos((y+6)*pi/180))/(cos((y+6)*pi/180)*cos((y+6)*pi/180))+sin(cos((x+7)*pi/180))*sin(cos((x+7)*pi/180))/(cos((x+7)*pi/180)*cos((x+7)*pi/180))*sin(cos((y+7)*pi/180))*sin(cos((y+7)*pi/180))/(cos((y+7)*pi/180)*cos((y+7)*pi/180))+sin(cos((x+8)*pi/180))*sin(cos((x+8)*pi/180))/(cos((x+8)*pi/180)*cos((x+8)*pi/180))*sin(cos((y+8)*pi/180))*sin(cos((y+8)*pi/180))/(cos((y+8)*pi/180)*cos((y+8)*pi/180)) +sin(cos((x+9)*pi/180))*sin(cos((x+9)*pi/180))/(cos((x+9)*pi/180)*cos((x+9)*pi/180))*sin(cos((y+9)*pi/180))*sin(cos((y+9)*pi/180))/(cos((y+9)*pi/180)*cos((y+9)*pi/180))+sin(cos((x+10)*pi/180))*sin(cos((x+10)*pi/180))/(cos((x+10)*pi/180)*cos((x+10)*pi/180))*sin(cos((y+10)*pi/180))*sin(cos((y+10)*pi/180))/(cos((y+10)*pi/180)*cos((y+10)*pi/180))+sin(cos((x+11)*pi/180))*sin(cos((x+11)*pi/180))/(cos((x+11)*pi/180)*cos((x+11)*pi/180))*sin(cos((y+11)*pi/180))*sin(cos((y+11)*pi/180))/(cos((y+11)*pi/180)*cos((y+11)*pi/180))+sin(cos((x+12)*pi/180))*sin(cos((x+12)*pi/180))/(cos((x+12)*pi/180)*cos((x+12)*pi/180))*sin(cos((y+12)*pi/180))*sin(cos((y+12)*pi/180))/(cos((y+12)*pi/180)*cos((y+12)*pi/180))+sin(cos((x+13)*pi/180))*sin(cos((x+13)*pi/180))/(cos((x+13)*pi/180)*cos((x+13)*pi/180))*sin(cos((y+13)*pi/180))*sin(cos((y+13)*pi/180))/(cos((y+13)*pi/180)*cos((y+13)*pi/180)) +sin(cos((x+14)*pi/180))*sin(cos((x+14)*pi/180))/(cos((x+14)*pi/180)*cos((x+14)*pi/180))*sin(cos((y+14)*pi/180))*sin(cos((y+14)*pi/180))/(cos((y+14)*pi/180)*cos((y+14)*pi/180))+sin(cos((x+15)*pi/180))*sin(cos((x+15)*pi/180))/(cos((x+15)*pi/180)*cos((x+15)*pi/180))*sin(cos((y+15)*pi/180))*sin(cos((y+15)*pi/180))/(cos((y+15)*pi/180)*cos((y+15)*pi/180))+sin(cos((x+16)*pi/180))*sin(cos((x+16)*pi/180))/(cos((x+16)*pi/180)*cos((x+16)*pi/180))*sin(cos((y+16)*pi/180))*sin(cos((y+16)*pi/180))/(cos((y+16)*pi/180)*cos((y+16)*pi/180))+sin(cos((x+17)*pi/180))*sin(cos((x+17)*pi/180))/(cos((x+17)*pi/180)*cos((x+17)*pi/180))*sin(cos((y+17)*pi/180))*sin(cos((y+17)*pi/180))/(cos((y+17)*pi/180)*cos((y+17)*pi/180))+sin(cos((x+18)*pi/180))*sin(cos((x+18)*pi/180))/(cos((x+18)*pi/180)*cos((x+18)*pi/180))*sin(cos((y+18)*pi/180))*sin(cos((y+18)*pi/180))/(cos((y+18)*pi/180)*cos((y+18)*pi/180))+sin(cos((x+19)*pi/180))*sin(cos((x+19)*pi/180))/(cos((x+19)*pi/180)*cos((x+19)*pi/180))*sin(cos((y+19)*pi/180))*sin(cos((y+19)*pi/180))/(cos((y+19)*pi/180)*cos((y+19)*pi/180)) title 'SHG for N=20'        
gnuplot> splot sin(cos(x*pi/180))*sin(cos(x*pi/180))/(cos(x*pi/180)*cos(x*pi/180))*sin(cos(y*pi/180))*sin(cos(y*pi/180))/(cos(y*pi/180)*cos(y*pi/180))+sin(cos((x+1)*pi/180))*sin(cos((x+1)*pi/180))/(cos((x+1)*pi/180)*cos((x+1)*pi/180))*sin(cos((y+1)*pi/180))*sin(cos((y+1)*pi/180))/(cos((y+1)*pi/180)*cos((y+1)*pi/180))+sin(cos((x+2)*pi/180))*sin(cos((x+2)*pi/180))/(cos((x+2)*pi/180)*cos((x+2)*pi/180))*sin(cos((y+2)*pi/180))*sin(cos((y+2)*pi/180))/(cos((y+2)*pi/180)*cos((y+2)*pi/180))+sin(cos((x+3)*pi/180))*sin(cos((x+3)*pi/180))/(cos((x+3)*pi/180)*cos((x+3)*pi/180))*sin(cos((y+3)*pi/180))*sin(cos((y+3)*pi/180))/(cos((y+3)*pi/180)*cos((y+3)*pi/180))+sin(cos((x+4)*pi/180))*sin(cos((x+4)*pi/180))/(cos((x+4)*pi/180)*cos((x+4)*pi/180))*sin(cos((y+4)*pi/180))*sin(cos((y+4)*pi/180))/(cos((y+4)*pi/180)*cos((y+4)*pi/180)) +sin(cos((x+5)*pi/180))*sin(cos((x+5)*pi/180))/(cos((x+5)*pi/180)*cos((x+5)*pi/180))*sin(cos((y+5)*pi/180))*sin(cos((y+5)*pi/180))/(cos((y+5)*pi/180)*cos((y+5)*pi/180))+sin(cos((x+6)*pi/180))*sin(cos((x+6)*pi/180))/(cos((x+6)*pi/180)*cos((x+6)*pi/180))*sin(cos((y+6)*pi/180))*sin(cos((y+6)*pi/180))/(cos((y+6)*pi/180)*cos((y+6)*pi/180))+sin(cos((x+7)*pi/180))*sin(cos((x+7)*pi/180))/(cos((x+7)*pi/180)*cos((x+7)*pi/180))*sin(cos((y+7)*pi/180))*sin(cos((y+7)*pi/180))/(cos((y+7)*pi/180)*cos((y+7)*pi/180))+sin(cos((x+8)*pi/180))*sin(cos((x+8)*pi/180))/(cos((x+8)*pi/180)*cos((x+8)*pi/180))*sin(cos((y+8)*pi/180))*sin(cos((y+8)*pi/180))/(cos((y+8)*pi/180)*cos((y+8)*pi/180)) +sin(cos((x+9)*pi/180))*sin(cos((x+9)*pi/180))/(cos((x+9)*pi/180)*cos((x+9)*pi/180))*sin(cos((y+9)*pi/180))*sin(cos((y+9)*pi/180))/(cos((y+9)*pi/180)*cos((y+9)*pi/180))+sin(cos((x+10)*pi/180))*sin(cos((x+10)*pi/180))/(cos((x+10)*pi/180)*cos((x+10)*pi/180))*sin(cos((y+10)*pi/180))*sin(cos((y+10)*pi/180))/(cos((y+10)*pi/180)*cos((y+10)*pi/180))+sin(cos((x+11)*pi/180))*sin(cos((x+11)*pi/180))/(cos((x+11)*pi/180)*cos((x+11)*pi/180))*sin(cos((y+11)*pi/180))*sin(cos((y+11)*pi/180))/(cos((y+11)*pi/180)*cos((y+11)*pi/180))+sin(cos((x+12)*pi/180))*sin(cos((x+12)*pi/180))/(cos((x+12)*pi/180)*cos((x+12)*pi/180))*sin(cos((y+12)*pi/180))*sin(cos((y+12)*pi/180))/(cos((y+12)*pi/180)*cos((y+12)*pi/180))+sin(cos((x+13)*pi/180))*sin(cos((x+13)*pi/180))/(cos((x+13)*pi/180)*cos((x+13)*pi/180))*sin(cos((y+13)*pi/180))*sin(cos((y+13)*pi/180))/(cos((y+13)*pi/180)*cos((y+13)*pi/180)) +sin(cos((x+14)*pi/180))*sin(cos((x+14)*pi/180))/(cos((x+14)*pi/180)*cos((x+14)*pi/180))*sin(cos((y+14)*pi/180))*sin(cos((y+14)*pi/180))/(cos((y+14)*pi/180)*cos((y+14)*pi/180))+sin(cos((x+15)*pi/180))*sin(cos((x+15)*pi/180))/(cos((x+15)*pi/180)*cos((x+15)*pi/180))*sin(cos((y+15)*pi/180))*sin(cos((y+15)*pi/180))/(cos((y+15)*pi/180)*cos((y+15)*pi/180))+sin(cos((x+16)*pi/180))*sin(cos((x+16)*pi/180))/(cos((x+16)*pi/180)*cos((x+16)*pi/180))*sin(cos((y+16)*pi/180))*sin(cos((y+16)*pi/180))/(cos((y+16)*pi/180)*cos((y+16)*pi/180))+sin(cos((x+17)*pi/180))*sin(cos((x+17)*pi/180))/(cos((x+17)*pi/180)*cos((x+17)*pi/180))*sin(cos((y+17)*pi/180))*sin(cos((y+17)*pi/180))/(cos((y+17)*pi/180)*cos((y+17)*pi/180))+sin(cos((x+18)*pi/180))*sin(cos((x+18)*pi/180))/(cos((x+18)*pi/180)*cos((x+18)*pi/180))*sin(cos((y+18)*pi/180))*sin(cos((y+18)*pi/180))/(cos((y+18)*pi/180)*cos((y+18)*pi/180))+sin(cos((x+19)*pi/180))*sin(cos((x+19)*pi/180))/(cos((x+19)*pi/180)*cos((x+19)*pi/180))*sin(cos((y+19)*pi/180))*sin(cos((y+19)*pi/180))/(cos((y+19)*pi/180)*cos((y+19)*pi/180))+sin(cos((x+20)*pi/180))*sin(cos((x+20)*pi/180))/(cos((x+20)*pi/180)*cos((x+20)*pi/180))*sin(cos((y+20)*pi/180))*sin(cos((y+20)*pi/180))/(cos((y+20)*pi/180)*cos((y+20)*pi/180)) title 'SHG for N=21' 

Second Harmonic Generation

gnuplot> set xrange [-180:180]
gnuplot> set yrange [-180:180]
gnuplot> set pm3d
gnuplot> set hidden3d 
gnuplot> set title 'SHG'
gnuplot> splot sin(cos(x*pi/180))*sin(cos(x*pi/180))/(cos(x*pi/180)*cos(x*pi/180))*sin(cos(y*pi/180))*sin(cos(y*pi/180))/(cos(y*pi/180)*cos(y*pi/180)) title 'N=1

'

Angle of Idler in Spontaneous Parametric Down Conversion

gnuplot> set xrange [-180:180]
gnuplot> set yrange [-180:180]
gnuplot> splot asin(sin(x*pi/180)*(y-400)/400)*180/pi 

FREE FALL


## A function, free fall, that takes in h (in metres), and
## returns the final velocity of the ball at the
## time step just before it touches the ground and makes
## a plot of velocity versus time.
function freefall(h)
## constants and initializations
v=0; ## initial velocity(at rest) [m/sec]
y=h; ## initial altitude of the ball from the ground [m]
t=0; ## initial time[sec ]
B1=0.05; ## Coeff of the term prop to v [kg/sec]
B2=6E-4; ## Coeff of the term prop to v^2 [kg/m]
m=0.25; ## Mass of the ball [kg]
g=9.8; ## Gravitational acceleration [m/sec^2]
dt=0.1; ## Time increment [sec]
n=1; ## Initialize the loop index[dimensionless]
## Run the loop or iterate until the vertical component is
## smaller than zero.
while (y(n)>0);
## decrease the altitude in each step
y=[y;y(n)+v(n)*dt]; ## and accumlate the results in the array
## of the vertical displacement.
## increase the time in each step
t=[t;t(n)+dt]; ## and accumlate the results in the array
## of the time for the time axis in the plot.
## increase the vertical velocity
v=[v;v(n)-dt*(g-B1*v(n)/m+B2*v(n)^2/m)]; ## in each step and
## accumlate the results
## in the array of vertical velocity.
n++; ## increase n by 1 in each step
endwhile
vfinal=-v(n) ## magnitude of the final velocity of the ball at the
## time step just before it touches the ground
plot(t,-v,';drag;') ## plot of the magnitude of the velocity versus time.
endfunction

INTEGRAL (TRAPEZOIDAL RULE)


## This m-file computes the integral of a function f over the
## equally spaced points x using the trapezoidal rule.
function I = trapezoid(x,f)
h = x(2)-x(1);
N = length(f);
I = 0;
for n=2:N-1,
I = I + h*f(n);
endfor
I = I + h/2*(f(1) + f(N));
endfunction

TOTAL ISING ENERGY


function E=Ising_Energy(lattice,i,j,rneigh)
H=2; ## applied magnetic field times permeability[joule]
E=0; ## initial configuration energy[joule]
J0=1 ## coupling constant of function J, normaly in joule.
alpha=2 ## exponential constant of function J.
L=length(lattice); ## number of particles or spins in lattice
for r=1:rneigh
NN1=mod(j+r,L); NN1+=(NN1==0)*L;
NN2=mod(j-r,L); NN2+=(NN2==0)*L;
NN3=mod(i+r,L); NN3+=(NN3==0)*L;
NN4=mod(i-r,L); NN4+=(NN4==0)*L;
J=coupling(rneigh,J0,alpha);
dE=J*(lattice(i,NN1)+lattice(i,NN2)+
lattice(NN3,j)+lattice(NN4,j))-H*(sum((rand(L)>0.5)*2-1));
E+=-1*dE
endfor
endfunction

NONUNFORM SOUND WAVE GENERATOR


## a script for a sound wave source made up of two different masses ##seperated in the middle. Let the rigth part is lighter than the
##left one that ensures that the ligther is faster, and so has
##greater r since r=cdt/dx which is dimensionless(ratio).
dx=1e-2; ## Increment in the horizontal distance (m)
c1=300; ## speed of the left one (m/s)
c2=500; ## Speed of the rigth one (m/s)
dt=dx/max(c1,c2); ## Increment in time (s). The max c, the min dt.
r1=c1*dt/dx; ## Ratio r for the left one
r2=c2*dt/dx; ## Ration r for the right one
x=-1:dx:1; ## same as the string_fixed
l=length(x);
x0=0.5;
k=1e2;
## Inital profile(return to guassuian distibution)
y=initial_profile(x,x0,k);
axis([-1.05,1.05,-1.1,1.1])
plot(x,y,'r;;',[0 0],[-1.1 1.1],'r-;;')
pause
## Boundary conditions while t=dtn=0
ynow=y;
yprev=y;
N=100;
for n=1:N
ynext=propagate_two_parts(ynow,yprev,r1,r2); ## calling the prepared
## function
plot(x,ynext,';;',[0 0],[-1.1 1.1],'r-;;')
axis([-1.05,1.05,-1.1,1.1])
pause(0);
yprev=ynow;
ynow=ynext;
endfor

STIFFNESS OF BENDING STRINGS


## A 'realistic', 'non-elastic' string, which responses to any
## bending and has stifness. This script takes in the previous
## and the present profiles and iterates to find
## the profile in the next time step. The ratio 'r' is not 1
## like in the 'non-realistic' string since the speed of the wave
## always less than the speed of the string it should be less than 1
## for best and most stable solution
##constants
dx=1e-2 ## Spatial increment (m)
L=2 ## Length of the string (m)
M=L/dx ## Dimensionless partition
r=0.25 ## Famous dimensionless ratio
E=1e-4 ## Dimensionless stiffnes
x=-1:dx:1;
l=length(x);
x0=0.5;
k=1e2;
## Set up the initial profile
y=initial_profile(x,x0,k);
plot(x,y)
pause
## Impose the time boundary condition
ynow=y;
yprev=y;
ynow(1)=ynow(2)=0
Nsteps=2000;
for n=3:Nsteps
ynow(Nsteps-1)=ynow(Nsteps-2)=0;
ynext=propagate_stiff(ynow,yprev,r);
plot(x,ynext,';;')
axis([-1.05,1.05,-1.1,1.1])
pause(0);
endfor

TRAVELING WAVE ON A NONREALISTIC STRING


## Script that simulates a traveling wave on an nonrealistic string.
##The one end of the string is driven sinusoidally
## while the other end is kept fixed.
dx=1e-2; ## increment in horizontal displacement(m)
c=300; ## speed of the wave (m/s)
dt=dx/c; ## Increment in time (s)
r=c*dt/dx; ## Ratio of speed of wave to
## the speed of the string
## Initial profile
x=-1:dx:1;
y=zeros(size(x)); ## Change the gaussian distribution
## to get zero vertical
##displacements for all parts of
##the string.(purely flat initially)
plot(x,y)
pause
##For the sinusoidally driven end, given constans are:
A=0.5; ##amplitude(m)
rate=100; ## inverse of the period (cycles/s)
omega=2*pi*rate ; ## angular frequency(Radians/s)
## Time boundary conditions. Get previous and
## present displacement as initial profile
ynow=y;
yprev=y;
Nsteps=2000;
## all steps are the same for the loop as the string_fixed
## apart from the function ynext
for n=1:Nsteps
ynext=propagate_driven(ynow,yprev,r,omega,A,n,dt);
plot(x,ynext,';;')
axis([-1.05,1.05,-1.1,1.1])
pause(0);
yprev=ynow;
ynow=ynext;
endfor

BIO-SAVART LAW


/*A C code that returns the magnetic field at a mesh of points
*on the xy-plane (xp,yp) for a straight wire segment through which
*current passes. The total magnetic field used in this code
*has the form (cos(theta1)-cos(theha2))/(yp*r_mag)
*ignoring the constants permeability, pi and current I and
*taking them together as 1. Thus no need for trapezoid since
*the integral of Bio-Savart can be tractable!!!
*r_mag represents the distance between the dl segment and
*the point (xp,yp) where we will find the B*/
#include <stdio.h>
#include <math.h>
/*Global constant*/
#define L 3.0 /* length of the wire*/
double r_mag(double x, double xp, double yp) {
double rx;
rx=xp-x; /* x is any point on xprime*/
return sqrt(rx*rx+yp*yp);
}
main()
{ int i;
double xp,yp;
double B; /*the z component of mag-field that is total field.*/
FILE *fid;
fid=fopen("mag-field","w");
/*loop over xp and yp values*/
for (yp=-2.0;yp<=2.0;yp=yp+0.1)
for (xp=-1.5;xp<=1.5;xp=xp+0.1)
{ /* Form the function B*/
B=2*(xp-L/2)/(yp*r_mag(L/2,xp,yp));
fprintf(fid,"%10.6lf %10.6lf %10.6lf\n",xp,yp,B);
}
fclose(fid);
}

TOTAL DISPLACEMENT OF LOOPLESS RANDOM WALK


## A for good and evil as well as loopless
## total displacement function of random walk
## that takes in the step numbers as an
## argument and returns the displacement
## of the rnd walker.Note that the probabilities
## to turn rigth and left are equal but the left
## step is twice the rigth step.
## Usage : rw_uneven(N)
function rw_uneven(N)
rn=rand(N,1); ## If the element of a random vector
r=-(rn<0.5); ## is smaller than 0.5, then return -1,
li=find(r == 0); ## Else return 2(equal probability).
r(li)=2;
x=sum(r) ## total displacement
endfunction

MEAN SQUARED DISPLACEMENT-1D RANDOM WALK


## Routine random walk(rw) for simulating
## one-dimensional random walk
## and calculating the mean squared displacement.
x2ave=[]; ## Initial array of the squared
## displacement at initial time(meters).
Nrw=3000; ## Number of walkers (Dimensionless).
beg=1000; ## Minimum step number ( // ).
inc=1000; ## Step number increment ( // ).
Nstepsmax=10000; ## Maximum step number ( // ).
dt=1; ## Time increment(second).
##Entering the nested loops.
for Nsteps=beg:inc:Nstepsmax; ## The loop through the number
Nsteps ## of steps each in walk.
x2=0; ## Initial squared location of all of the
## walkers at all steps.
for m=1:Nrw ## The loop through the desired
## number of walkers.
r=rand_disc_rev(Nsteps,0.5); ## Calling the rw generator
## function which will give column vector
## of each elements are either 1 or -1.
x2+=sum(r)^2; ## Total displacement squared whose
## elements stems from the vector r.
endfor
x2ave=[x2ave;x2/Nrw]; ## Accumulation of the squared displacements
endfor ## at the step Nsteps into the array x2ave(Nsteps).
## Exiting the nested loops. Continuing by constructing polynomials.
hold off
plot(dt*[beg:inc:Nstepsmax],x2ave,'b*;RawData;') ## Plot of the raw data
## of mean squared displament vs. step
## number or here the time.
[p,s]=polyfit(dt*[beg:inc:Nstepsmax]',x2ave,1); ## Calculate the coefficients
## p and the quality of measure s
## of the 'imaginary' polinomial, x2ave(Nsteps).
pval=polyval(p,dt*[beg:inc:Nstepsmax]); ## The values of the fit polynomial
## in the time interval dt*[beg,Nstepsmax].
hold on
plot(dt*[beg:inc:Nstepsmax],pval,'r-;Fit;') ## Plot of fit data of x2ave vs time.
hold off

RANDOM WALK GENERATOR


## This function
## provides us with very simple
## random walk generator compared
## to the function rand_disc_loop
## which is time-consuming, in operation,
## long, and unuseful for other than the
## coin toss simulatons
function xy=rand_disc(N); ## RW generator.
r=floor(rand(N,1)*4); ## Random coulumn vector of N
## integer elements
## multiplied by 4 to widen the
## interval from [0,1] to [0,3].
x=y=zeros(size(r)); ## The coulumn vectors
## of N elements with all elements zero.
x(find(r==0)) = 1; ## The elements of x,
## the vector function
## which takes in the row #
## of the zero elements of
## the vector r as an argument,
## is assigned to 1.
x(find(r==1)) =-1; ## The elements of x,
## the vector function
## which takes in the row # of the
## elements of 1 of the vector r
## as an argument,
## is assigned to -1.
y(find(r==2)) = 1; ## The elements of y,
## the vector function
## which takes in the row # of the
## elements of 2 of the vector r
## as an argument,
## is assigned to 1.
y(find(r==3)) =-1; ## The elements of y,
## the vector function
## which takes in the row # of the
## elements of 3 of the vector r
## as an argument,
## is assigned to -1.
xy=[x y]; ## The resulting x by y
## matrix of a random walk
## with elements 1 and -1.
endfunction

PROPAGATION OF NON-UNIFORM STRING


## The string is made up of two different mass density strings
## seperated in the middle. Choose the right one to be ligther
## than the left. So, since velocity is inversely propotional
## to the mass density, than the rigth one is faster also.
## r1=c1*dt/dx and r2=c1*dt/dx
function ynext=propagate_two_parts(ynow,yprev,r1,r2)
## initiallization for, take y as previos one
ynext=ynow;
## we have two parts, hence two loops are needed
## the first part is between x=0 and
## let y(now)/2 for better visiulation.
## floor(x) returns the largest integer not greater than x
## length(x) determines the number of column
## or rows in matrix or vector
##I started the 'for' from x=0 but didn't work
##then started from x=1
## but in this case the
## string_two_parts didn't work
## let x0=i0=2
for i=2:floor(length(ynow)/2)
ynext(i) = 2*(1-r1^2)*ynow(i)-yprev(i)+r1^2*(ynow(i+1)+ynow(i-1));
endfor
## the second part is between x0=(L1+L2) /2 and L1+L2 where
## L1 and L2 are string lenghts respectively
## I started the 2nd for from N2steps
## and finished it by 'ynow' but it
## didn't work. So i tried the following
for i=floor(length(ynow)/2)+1:floor(length(ynow))-1
ynext(i) = 2*(1-r2^2)*ynow(i)-yprev(i)+r2^2*(ynow(i+1)+ynow(i-1));
endfor
endfunction

REALISTIC AND NON-REALISTIC STRINGS


## A 'realistic', 'non-elastic' string, which responses to any
## bending and has stifness. This script takes in the previous
## and the present profiles and iterates to find
## the profile in the next time step. The ratio 'r' is not 1
## like in the 'non-realistic' string since the speed of the wave
## always less than the speed of the string it should be less than 1
## for best and most stable solution
##constants
dx=1e-2 ## Spatial increment (m)
L=2 ## Length of the string (m)
M=L/dx ## Dimensionless partition
E=1e-4 ## Dimensionless stiffnes
function ynext=propagate_stiff(ynow,yprev,r)
## Quick and dirty way to fix boundary conditions -- for each step
## they are the same as the previous step.
ynext=ynow;
ynow(1)=ynow(2)=0;
##Entering the loop
for i=3:length(ynow)-1
## boundary condition
ynow(length(ynow)-1)=ynow(length(ynow)-2)=0;
## Divide the ynext with many terms into three parts for easiness
ynext(i)=(2−(2*r^2)−(6*E*(r^2)*(M^2)))*ynow(i)−yprev(i);
ynext(i)=ynext(i)+(r^2)*(1+4*E*(M^2))*(ynow(i+1)+ynow(i−1));
ynext(i)=ynext(i)-E*(r^2)*(M^2)*(ynow(i+2)+ynow(i−2));
endfor
endfunction

SINUSOIDALY DRIVEN STRING


## Takes in the previous and the present profiles of the string and
## iterates to find the profile in the next time step.
function ynext=propagate_driven(ynow,yprev,r,omega,A,n,dt)
## initialization, boundary condition
ynext=ynow;
## takes in length(ynow) as an argument
## since one end of the string is driven sinuosidaly
## omega is the angular frequency of the driven force
## n is the time index, n*dt is t(n)
ynext(length(ynow))=A*sin(omega*n*dt);
## entering the loop
for i=2:length(ynow)-1
ynext(i) = 2*(1-r^2)*ynow(i)-yprev(i)+r^2*(ynow(i+1)+ynow(i-1));
endfor
endfunction

PRIME OR NOT


## A funtion that takes a number
## in as an argument and returns
## the output p=1 if it is prime
## otherwise returns p=0.
function prime_or_not(num)
p=1; ## assume the number is prime.
div=2; ## the smallest divisor.
## since 2*num/2> num. Forbidden
Nsteps=num/2 ; ## Interval of division: the half
## of the num end itself
for i= div:Nsteps ## to divide al the numbers less than
## num/2.
if(rem(num,i)==0) ## any prime num cannot be
## divided with an integer.
p=0; ## otherwise exit the loop.
break;
endif
endfor
if(p == 1) ## prints the result
p=1
else
p=0
endif
endfunction

FOURTH ORDER RUNGE-KUTTA METHOD- MOTION OF A SPHERICAL MASS WITH AIR RESISTANCE


## Motion of a spherical mass with air resistance
## Fourth order Runge-Kutta Method
## Very Important!!! The positive velocity direction
## is the direction of the gravitational acceleration
m=1E-2; ## mass of the object(kg/m)
g=9.8; ## Acceleration due to the gravity (m/sec^2)
v0=0; ## initial velocity of the object(m/sec)
k=1E-4; ## Air drag coefficient(kg/m)
t=[]; ## Empty time vector(sec)
t(1)=0; ## Released time(sec)
dt=1E-1; ## Increment in time(sec)
v=[]; ## Empty velocity vector(m/s)
v(1)=v0; ## Initial velocity
v_nodrag=[]; ## Velocity by ignoring air drag
## Analytic solution by zeroth order approximation.
v_nodrag(1)=0; ## Velocity vector with no drag.
n=1; ## Initialize the loop index
f0=[];
f1=[]; ##The 1st order derivatives by
f2=[]; ##by Runge-Kutta. Eqn 5.33 pg217
f3=[];
## Run the loop until the time reaches the value 10sec.
while (t(n)<=10);
f0(n)=g-(k/m)*v(n)*v(n);
v_f0(n)=v(n)+(dt/2)*f0(n);
f1(n)=g-(k/m)*v_f0(n)*v_f0(n);
v_f1(n)=v(n)+(dt/2)*f1(n);
f2(n)=g-(k/m)*v_f1(n)*v_f1(n);
v_f2(n)=v(n)+dt*f2(n);
f3(n)=g-(k/m)*v_f2(n)*v_f2(n);
v(n+1)=v(n)+(dt/6)*(f0(n)+2*f1(n)+2*f2(n)+f3(n));
t(n+1)=t(n)+dt;
v_nodrag(n+1)=v_nodrag(n)+g*dt; ##analytic solution
n++;
endwhile
plot(t,v,'r-',t,v_nodrag,'b-');
title('Velocity vs Time');
xlabel('time(sec)');
ylabel('velocity(m/sec)');
legend('v(Runge-Kutta)','v(no drag)')
axis([0,13]);
save -text RUNGEKUTTA.dat
print('-dpsc','RUNGEKUTTA.ps ')

JACOBI AND GAUSS-SIEDEL SCHEMES


##The formulas and the iteration procedure
##below are taken from the Devies Book
##Eqn 5.89
##Constants and Initializations
dx=1E-3; ##Increment in x.
c1=1-2.5*dx; ## c1,c2,c3,c4 are constants in Eqn.5.89
c2=1+2.5*dx;
c3=-10*dx*dx;
c4=1/(2+c3);
## Set up the initial profile
x=0.0:0.01:1.0;
yold=0:100;
hold on
plot(x,yold);
xlabel('x');
ylabel('y');
##Impose the boundary conditions
Nsteps=1000;
yoldJac=yold; ## old y values for Jacobi Iteration Scheme
yoldGauSied=yold; ## old y values for Gauss-Siedel Scheme
ynewJac=yoldJac; ## new y values for Jacobi Iteration Scheme
ynewGauSied=yoldGauSied; ## new y values for Gauss-Siedel Scheme
for n=1:Nsteps
for i=2:100
ynewJac(i)=c4*(c1*yoldJac(i+1)+c2*yoldJac(i-1)+c3*x(i));
ynewGauSied(i)=c4*(c1*yoldGauSied(i+1)+c2*ynewGauSied(i-1)+c3*x(i));
endfor
yoldJac=ynewJac;
yoldGauSied=ynewGauSied;
endfor
plot(x,ynewJac,'g*',x,ynewGauSied,'m@');
title('y vs x');
xlabel('x');
ylabel('y');
legend('y(initial)','Jacobi','Gauss-Siedel');
axis([0,1.5]);
hold off
save -text JACOBIGAUSSSIEDEL.dat
print ('-dpsc','JACOBIGAUSSSIEDEL.ps')

SIMPLE AND MODIFIED EULER METHOD


##Usage:Call Octave from terminal
##and then call EulerMethodUmitAlkus.m
##from octave and finally
##press enter. That's all.
##Simple Euler Method
##Constants and initializations
x=[]; ## initial empty vector for x
y=[]; ## initial empty vector for y
x(1)=1; ## initial value of x
y(1)=1; ## initial value of y
h=1E-3; ## increment in x
dery=[]; ## 1st derivative of y wrt x
dery(1)=0;## 1st entry of dery
n=1; ## inital loop index for while
## enter the while loop for the interval x=[1,2]
while (x(n)<=2)
x(n+1)=x(n)+h;
dery(n+1)=x(n)*x(n)-2*y(n)/x(n); ##given
y(n+1)=y(n)+h*dery(n+1); ##Euler method
n++;
endwhile
##exit from the 1st while loop
##Modified Euler Method
##Constant and initializations
x(1)=1; ## beginnig of the interval [1,2]
ymod(1)=1; ## inital value for modified y.
ymid=[]; ## empty vector function evaluated at x midpoint
xmid=[]; ## empty vector func. of midpoints of the interval h in x-axis.
derymod=[]; ## modified derivatives of ymod wrt x.
derymid=[]; ## derivative of ymid wrt xmid.
derymid(1)=0; ##1st entry of der. of y wrt x at midpoints.
##Enter the 2nd while loop
n=1;
while (x(n)<=2) ## x(n)'s are the beginning values of the interval h.
xmid(n)=x(n)+h/2;
derymod(n)=x(n)*x(n)-2*ymod(n)/x(n); ##given equation
ymid(n)=ymod(n)+derymod(n)*h/2;
derymid(n+1)=xmid(n)*xmid(n)-2*ymid(n)/xmid(n);
ymod(n+1)=ymod(n)+derymid(n)*h; ##modified Euler Method
x(n+1)=x(n)+h;
n++;
endwhile
## Plot for Simple Euler Method
subplot(2,1,1)
hold on
plot(x,dery,'c-');
plot(x,ymod,'r-');
legend('dy/dx','y')
title('Simple Euler Method');
xlabel('x');
hold off
##Plot for Modified Euler Method
subplot(2,1,2)
hold on
plot(x,derymid,'m-');
plot(x,ymod,'g-');
title('Modified Euler Method');
xlabel('x');
legend('dy/dx(mid)','y-mod');
hold off
print('-dpsc','SIMMODEULER.ps')
save -text SIMMODEULER.dat

NEWTON-RAPSON METHOD-8th degree Legendre polynomial


## Newton-Rapson Method to the smallest non negative root
## of the 8th degree Legendre Polynomial
## P8(x)=(1/128)(6435x^8-12012x^6+6930x^4-1260x^2+35)
## where -1<=x<=1.
## for the smallest non negative root, we can ignore
## all the terms except the last two by truncated
## the function to be zero and find
## x=0.167 as the initial smallest non negative
## root.
##Constants and initializations
x=[]; ## Empty array for the iterated x roots
x(1)=0.16700000; ## Initial guess to begin the iteration for the
## smallest non-negative root.
L8=[]; ## Empty array for the Legendre polynomial
L8p=[]; ## Empty array for the derivative of the Legendre polynomial
for i=1:100
##The value of the function at x
L8(i)=(1/128)*(6435*x(i)^8-12012*x(i)^6+6930*x(i)^4-1260*x(i)^2+35);
##The value of the derivative of the function at x
L8p(i)=(1/128)*(6435*8*x(i)^7-12012*6*x(i)^5+6930*4*x(i)^3-1260*2*x(i));
x(i+1)=x(i)-L8(i)/L8p(i); ## the iteration
endfor
## For plot let's define a new variable t instead of x and P8(t) instead of L8
t=-1:0.01:1;
P8=(1/128)*(6435*t.^8-12012*t.^6+6930*t.^4-1260*t.^2+35);
plot(t,P8);
title('8th degree Legendre polynomial vs x');
xlabel('x');
ylabel('P8(x)');
legend('P8(x)');
printf('The smallest non-negative root of the 8th Legendre polynomial is=%f',x(i))
save -text NEWRAPLEGENDRE.dat
print('-dpsc','NEWRAPLEGENDRE.ps');

MINIMUM OF GAMMA FUNCTION

## This string finds the minimum of the Gamma function
## integral(x^(alpha-1)*exp(-x))dx) (x=[0,infinity])
## in the interval 0<alpha<4.
dx=1E-2; ## Increment in x
x=0:dx:100; ## x array
Gamma=[]; ## Empty Gamma function
alpha=0.00:dx:4.00; ## The independent variable of Gamma function
for i=1:length(alpha)
## Call the trapezoid function to
## calculate all entries of Gamma function
## corresponding to entries of alpha
f=x.^(alpha(i)-1).*exp(-x); ## integrand of Gamma function
Gamma(i)=trapezoid(x,f);
endfor
plot(alpha,Gamma);
title('Gamma values vs alpha');
legend('Gamma(alpha)');
xlabel('alpha');
ylabel('Gamma');
## The following 'for loop' finds the minimum value of Gamma
## and the corresponding alpha value.
minimum=Gamma(1);
for n=1:length(Gamma)
if(minimum>Gamma(n))
minimum=Gamma(n);
n; ## The array index where the minimum of Gamma is.
m=alpha(n); ## The alpha value corresponding to the min of Gamma
endif
endfor;
printf('Minimum of the Gamma Function in the interval 0<alpha<4 is %f\n for the alpha=%f\n' ,minimum,m)
print('-dpsc',MINGAMMA.ps');
save -text MINGAMMA.dat

NEWTON-RAPSON METHOD FOR HEAT FLOW


##Constants and initializations
a=5.67E-8; ## Stefan-Boltzman constant[Watt/meter^2Kelvin^4]
e=0.8; ## Rod surface emissivity [Dimensionless]
h=20; ## Heat transfer coefficient of air flow [W/m^2-K]
Tinf=Ts=25; ## Temperature of air and the walls of the close[Celcius]
D=0.1; ## Diameter of the rod[meter]
I2R=100; ## Electric power dissipated in rod (Ohmic Heat)[W]
T=[]; ## Temperature of the rod[*C]
T(1)=25; ## Initial guess of the temperature of the rod[*C]
Q=[]; ## Heat function [W]
Qp=[]; ## First derivative of Q wrt T [W/C*].
for i=1:100
Q(i)=pi*D*(h*(T(i)-Tinf)+e*a*(T(i)^4-Ts^4))-I2R;
Qp(i)=pi*D*(h+4*e*a*T(i)^3);
T(i+1)=T(i)-Q(i)/Qp(i); ## Newton-Rapson Method
endfor
printf('The steady state temperature is %f\n',T(i+1))
save -text HeatFlowTemp.dat
## The plot
t=1:100; ##temperature
for n=1:100
H(n)=pi*D*(h*(t(n)-Tinf)+e*a*(t(n)^4-Ts^4))-I2R;
endfor
plot(t,H)
xlabel('T(Celcius)');
ylabel('Q(Watt)');
legend('Q(T)');
title('Heat flow vs Temperature')
print('-dpsc','HeatFlowTemp.ps')

FACTORIAL FOR KEPLER’S MOTION


##Constants and initializations
ecc=0.1; ## Eccentricity
M=24.851090; ## Mean anomaly(degrees)
Mr=M*pi/180; ## Mean anomaly (rads)
ser=[]; ## Series in E
jser=[]; ## Bessel function
for m=1:100
for n=0:100
jser(n+1)=((-1)^n)*((m*ecc/2)^(2*n+m))/(factorial(n)*factorial(m+n));
jser(n+1)+=jser(n+1);
endfor
ser(m)=(1/m)*jser(n+1)*sin(m*Mr);
ser(m)+=ser(m);
endfor
E=M+2*ser(m);
save -text FACTORIALKEPLER.dat
printf('We find E to be %f. This is the answer of Q1(c)\n',E);

BESSEL FOR KEPLER’S MOTION

##Constants and initializations
ecc=0.1; ## Eccentricity
M=24.851090; ## Mean anomaly(degrees)
Mr=M*pi/180; ## Mean anomaly (rads)
ser=[]; ## Series in E
for m=1:1000
ser(m)=(1/m)*besselj(1,m*ecc)*sin(m*Mr);
ser(m)+=ser(m);
endfor
E=M+2*ser(m);
save –text BESSELKEPLER.dat
printf('We find E to be %f. This is the answer of Q1(b)\n',E);

NEWTON-RAPSON METHOD FOR KEPLER’S ECCENTRICITY

##Constants and initializations
ecc=0.1; ## Eccentricity
M=24.851090; ## Mean anomaly(degrees)
Mr=M*pi/180; ## Mean anomaly (rads)
f=[]; ## The fundamental eqn we want to solve
## Kepler's Equation's zero function f(E)=E-esinE-M=0
fp=[]; ## 1st derivative of f(E). fp(E)=1+ecosE
Ed=[]; ## Eccentric anomaly(degrees)
Er=[]; ## Eccentric anomaly(rads)
Ed(1)=M+0.85*sign(sin(Mr))*ecc; ## The eccentric anomaly initial
## guess for Danby's method(degrees)
Er(1)=Ed(1)*pi/180; ## rads
for i=1:100
f(i)=Ed(i)-ecc*sin(Er(i))-M;
fp(i)=1+ecc*cos(Er(i));
Ed(i+1)=Ed(i)-f(i)/fp(i); ## The iteration for an updated eccentric anomaly E(i+1)
## based on a current value E(i) by Newton-Rapson Method
Er(i+1)=Ed(i+1)*pi/180;
endfor
save -text NEWTONRAPSON.dat
printf('We find E to be %f. This is the answer of Q1(a)\n',Ed(i+1));

NEWTON’S METHOD FOR MINIMUM


##Newton's Method to find
##the minimum of the function F(x)=(x-2)^4-9
##with the initial guess xmin=1.0
##Constants and initializations
xmin=[]; ##The empty array of x that minimizes the F(x)
xmin(1)=1.0; ##Initial value of the xmin
Fmin=[]; ##Minimum values of F(x)
x=0.0:0.1:4.0; ##Only for plotting purposes
F=[]; ##Our examined Function evaluated on x-space
Fp=[]; ##First derivative of F(x) wrt x
Fpp=[]; ##Second derivative o F(x) wrt x
NSteps=50; ##Step number of iteration
##Algorithm
for n=1:NSteps
Fmin(n)=(xmin(n)-2)^4-9;
Fp(n)=4*(xmin(n)-2)^3;
Fpp(n)=12*(xmin(n)-2)^2;
xmin(n+1)=xmin(n)-Fp(n)/Fpp(n);
Fmin(n+1)=(xmin(n+1)-2)^4-9;
endfor
printf("x*, at which F(x) is minimum, is %1.6f\n",xmin(n+1))
printf("Minimum of F(x) is %1.6f\n",Fmin(n+1))
F=(x-2).^4-9;
subplot(2,1,1)
plot(x,F)
title('Newton^,s Method-F(x) vs x');
xlabel('x');
ylabel('F(x)');
text(2,-7,'\downarrow')
text(1.7,-5.6,'(xmin,Fmin)=(2,-9)')
axis([-0,4,-10,10])
subplot(2,1,2)
plot(xmin,Fmin,'r-')
title('F(xmin) vs xmin')
xlabel('xmin')
ylabel('F(xmin)')
axis([0,4,-10,-7])
save -text NEWMETUmitAlkus.dat
text(1,-9.6,'\downarrow')
text(0.6,-9.4,'initial guess for xmin')
print('-dpsc', 'NEWMETUmitAlkus.ps')

NUMBER OF CLUSTER POINTS

#include <stdio.h>
#include <math.h>
main()
{
int i;
double points[99][2];
FILE *fid,*write;
fid=fopen("coords.txt","r");
for(i=0;i<100;i++)
points[i]={x,y,z};
fscanf(fid,"%1f%1f%1f",points[i]);
if(sqrt(x*x+y*y+z*z)==i)
write=fopen("out.txt","w");
fprintf(fid,"%d%d\n",i,points[i]);
fclose(write);
else
break;
fclose(fid);
}

INITIAL POSITIONS OF DIFFUSIVE PARTICLES

## A function takes in the number of diffusive particle
## as an argument and gets the initial positions in 2D
## as an outcome by inserting them in a square.
function p=initial_positions7(Nwalkers)
dt=3; ## to extend the interval of axes.
a=sqrt(Nwalkers); ## step size.
n=(a-1)/2; ## max range in the x and y axes.
if(rem(a,2)~=1) ## Nwalkers must be odd squared.
printf("The number you have entered is not an odd squared. Exiting...\n");
return
endif
n=(a-1)/2; ## max range in the x and y axes.
x=[-n:n]'; ## integer interval in the x axis.
p=zeros(a,2); ## initialize the positions of 'a' walkers.
p(:,1)=x; ## equate the column of p to x vector.
p=repmat(p,a,1); ## replicate positions of 'a' walkers to get a*a walkers.
for i=1:a:Nwalkers ## the loop through the number walkers.
m=(i+a-1)/a; ## index get the entries of x resp.
p(i:i+a-1,2)=x(m) ; ## y components are a-fold degenerate.
endfor
plot(p(:,1),p(:,2),'b*;;')
axis(dt*[-n, n,-n,n])
endfunction

FACTORIAL


## Function that calculates the factorial of a number
## Usage : f=factorial(n)
function f=factorial(n)
## Initialize the output
f=1;
## Check whether the input is correct
if ( (n<0) || (rem(n,1)~=0) )
printf("n cannot be a negative number. Exiting...\n");
return
endif
for num=1:n
f*=num;
endfor
endfunction

Simple Euler Method

##Usage:Call Octave from terminal
##and then call EulerMethodUmitAlkus.m
##from octave and finally
##press enter. That's all.
##Simple Euler Method
##Constants and initializations
x=[]; ## initial empty vector for x
y=[]; ## initial empty vector for y
x(1)=1; ## initial value of x
y(1)=1; ## initial value of y
h=1E-3; ## increment in x
dery=[]; ## 1st derivative of y wrt x
n=1; ## inital loop index for while
## enter the while loop for the interval x=[1,2]
while (x(n)<=2)
x(n+1)=x(n)+h;
dery(n+1)=x(n)*x(n)-2*y(n)/x(n); ##given
y(n+1)=y(n)+h*dery(n); ##Euler method
n++;
endwhile
##exit from the 1st while loop
##Modified Euler Method
##Constant and initializations
ymid=[]; ## empty vector function evaluated at x midpoint
xmid=[]; ## empty vector func. of midpoints in x
ymid(1)=1; ## inital value for ymid.
derymid=[]; ## derivative of y at midpoints
##Enter the 2nd while loop
n=1;
while (x(n)<=2)
xmid(n)=x(n)+h/2;
derymid(n)=xmid(n)*xmid(n)-2*ymid(n)/xmid(n);
ymid(n+1)=ymid(n)+derymid(n)*h/2;
n++;
endwhile
## Plot for Simple Euler Method
subplot(2,1,1)
hold on
plot(x,y,'r-');
title('Simple Euler Method');
xlabel('x');
plot(x,dery,'b-');
legend('y','dy/dx')
hold off
##Plot for Modified Euler Method
subplot(2,1,2)
plot(x,ymid,'g-');
title('Modified Euler Method');
xlabel('x');
ylabel('y-modified');
print('-dpsc','EulerMethodUmitAlkus.ps')
save -text EulerMethodUmitAlkus.dat

DIFFUSION-RANDOM WALK

## the function that takes in the initial positions
## of random walkers and the number of steps as
## an argument and returns the final positions
## in 2D. The figure of the walkers after the comple-
## tion of the walk also plotted.
function f=diffusion(Nsteps)
f=[];
Nwalkers=9;
for i=1:Nsteps ## The loop through the number
## of steps each in walk.
p=initial_positions(Nwalkers); ## Initial squared location of all of the
## walkers at all steps.
for m=1:Nwalkers ## The loop through the desired
## number of walkers.
r=rand_disc_rev(N); ## Calling the rw generator
## function which will give column vector
f=[f;p+r] ## of each elements are either 1 or -1.
## Total displacement whose
## elements stems from the vector r.
endfor
f
endfor
plot(f(:,1),f(:,2),'r*;Dif;')
endfunction

ISING COUPLING

## interaction energy for ising system between the spins seperated with distance r.
## J0 and alpha can be varied with respect to the ising system
function J=coupling(r,J0,alpha)
J=J0*exp((-r-1)/alpha);
Endfunction

Equation of Motion of Proton (Coulomb Potential)

## A function of the solution of the path
## equation of a proton under the inverse
## square atraction field of an electron
## that takes in the initial seperation distance
## or the position r0 wrt center of the electron
## as an argument and returns the total time
## ttotal it takes for it to reach within 1.0m
## of the electron and plots the velocity v
## vs time graph. usage ttotal=coulomb(r0).
function coulomb(r0)
## constants and initializations
k=9e+9; ## coulomb field constant [Nm^2/C]
q=1.6e-19; ## electronic charge [C]
mp=1.7e-27; ## proton mass [m]
dt=1e-4;; ## increment in time [sec]
r=r0; ## initial seperation[m]
t=0; ## initial time [s]
v0=0; ## initial velocity [m/s]
v=v0; ## 1st entry of v array [m/s]
n=1; ## initialization of loop index
## since we have already n=0 in argument r0.
while(r(n)>0.1);
dr=v(n)*dt; ## increment that is decrease in r
## since v(n) will be negative below.
r=[r;r(n)+dr]; ## decreases r in each step and
## accumulates the results in r array.
dv=-k*q*q*dt/(mp*r(n)*r(n)); ## increment in v to make v more negative.
v=[v;v(n)+dv]; ## increases the magnitude of negative v in each step and
## accumulates the results in v array.
t=[t;t(n)+dt]; ## increases the time in each step and
## accumulates the results in t array
## for both time axis and ttotal.
n++; ## increase n by 1 in each step
endwhile
ttotal=t(n) ## print the last entry of t array gives total time.
plot(t,v,';;')
xlabel('time(sec)'); ## Ola yuppi! I've learned eventually
ylabel('velocity(m/sec)'); ## to label the axes :) happy end!
endfunction

One Dimensional Harmonic Oscillator-Numerov Method

x=[];
h0=1;
M=4;
N=M+1;
x(1)=0;
x(N)=x(1)+h0*M;
x=x(1):h0:x(N)
A=zeros(N);
A(1,1)=-2*(5*(x(N)*h0)^2/12+1);
A(N,N)=-2*(5*(x(1)*h0)^2/12+1);
A(1,2)=1-(x(M)*h0)^2)/12;
A(N,M)=1-(x(2)*h0)^2)/12;
B=zeros(N);
B(1,1)=B(N,N)=-10*(h0^2)/6;
B(1,2)=B(N,M)=-(h0^2)/6;
for i=2:M
B(i,i)=-10*(h0^2)/6;
B(i,i-1)=B(i,i+1)=-(h0^2)/6;
A(i,i)=-2*(5*(x(N+1-i)*h0)^2/12+1);
A(i,i+1)=1-(x(N-i)*h0)^2)/12;
A(i,i-1)=1-(x(N+2-i)*h0)^2)/12;
end
A
B