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

Sunday, August 17, 2014

PHYSICS MACHINE

Physics Machine 
Ümit Alkuş

 Abstract
Physics machine is a software which does physics like a physicist. First, all the things human being has developed so far, for doing physics, will be available to this machine. Secondly, all the consistent theories, successful experiments, and published articles will be included into this machine in the form of traced and readable knowledge, in other words, this machine can read and understand these all. Finally, as the last target, this machine can observe the universe and physical events with the aim of creating theories and physical laws.

 METU, Physics Department, 06800, Ankara, Turkey
  Keywords: Artificial Intelligence, Machine Learning, Data Mining, Artificial Physicist

  Introduction
There are approximately millions of articles over physics, huge collection of very successful theories, and physics books. In the earth, no physicist could have attempted to read and understand these accumulations since it requires free infinite time, free space, infinite energy, and feeding from also other disciplinaries and fields. Especially, philosophy, science history and mathematics are needed in a very high understanding. Therefore, the first purpose of this thesis, which is to inject all the needs to do physics, appears here to tackle with these vital humanity problems. Hence, this research can be thought as an intention to meet all these impossibilities.

  Problems of Humanity
Human being has the age limit which is approximately 120; however intellectual people are aging faster than others. Of course, this handicaps the productive physicists from doing much more physics, especially in the very bright years of her/him . Therefore, the need for this physics machine seems to one obvious since it is immortal. On the other hand, this machine is not only older but also more productive than humanity. It consists, recombines and drives artificial physicists of all over the universe (no need for limited to the earth), physics knowledge, and physical productions. Then, this machine can operate like an artificial physical institution. Probably, this machine will be able to give us the key to develop memory transfer applicable to humanity. Thus, the subsequent mission of this physics machine will be the biology machine. So far in the history of humanity, unfortunately, no one, no state, no institution could have brought these things regarding physics together in the real world and this seems to be impossible for the future humanity because there is no attemption to solve this vital problem. Knowledge is infinite. No one person or even very powerful states, e.g. US, CCCP, or EU can own these infinite knowledge. For this, humanity needs a huge machine. Therefore, this physics machine can be considered to be an initiative huge machine. 
Humanity is like a balloon which is filled with untraced and unreadable knowledge gas. However, the developments especially in computer science and philosophy has the ability of giving the key to structure and invent this software machine. Some of the important fields in computer science which can be used in the construction of this physics machine are Artificial Intelligence, Machine Learning, and Data Mining. In the following lines, the researcher will analyze these fields so that they can be used to implement this machine.

  Artificial Intelligence
 Artificial intelligence is a type of intelligence which has been the subject of not only computer science but also philosophy. It can be regarded as the intelligence of nonliving organisms. According to the researcher, since the organic intelligence of human being, the institutions, the universities, or the states have no attempt to mine and gather physics knowledge to establish much more vital physical laws for superhumanity; the inorganic and nonliving intelligence, that is artificial intelligence is the best candidate to operate this machine.

 Machine Learning
Like humanity and its efforts to learn so far in the education systems, this physics machine is ready to learn and grow. Learning in computer is easier than that in humanity since to trace and store the knowledge, transmit information, sorting and memory transfer are faster and permanent than those in human brain. The developments in machine learning field are used to do all these processes so that this physics machine can operate the knowledge and practice to learn new things in physics.

  Data Mining
The most important thing to construct a theory in physics is to collect data as the result of intense physical observations. In addition, to extract the correct and usable data has vital importance. Data mining, then, will become the other important operator to construct the physics machine.

  Artificial Physicist
To gather human physicists together in order to construct advanced machines or theories is very difficult in this labor division system all over the world. Physics and physicist are divided into million parts in terms of field and this is the largest handicap to establish a large physical laws. This is because they cannot combine their knowledge to produce new things. Even after this automation of their knowledge and experiences, they leave with their previous knowledge and practices, so they do not learn new things. Thus, artificial physicist can solve this division problem in the construction of the physics machine.

Friday, May 16, 2014

LETTERS TO PROF. DR. WILLIAM WOOTTERS ON QUANTUM MIXED STATES

FIBER OPTIC EXPERIMENT II

Thank you Assoc. Prof. Dr. Hakan Altan

FIBER OPTIC EXPERIMENT I

Thank you Assoc. Prof. Dr. Hakan Altan

REVIEW OF QUANTUM MECHANICS II FINAL EXAM

Thank you Prof. Dr. Namık Kemal Pak

DECOMPOSITION OF DIBORANE (B2H6) AND ADSORPTION OF BH3 ON GE(100)

Thank you Prof. Dr. Şenay Katırcıoğlu

TOPICS IN QUANTUM INFORMATION THEORY EXAM SOLUTIONS

Thank you Assoc. Prof. Dr. Sadi Turgut

QUANTUM CRYPTOGRAPHY HOMEWORK SOLUTIONS

Thank you Prof. Dr. İbrahim Yurdahan Güler

Advanced Quantum Mechanics Homeworks Solutions

Thank you Prof. Dr. Cemal Yalabık

Fundamentals of Physics I-II ---My Lecture Videos- Have Fun!

Thank you Assoc. Prof. Dr. Hüseyin Oymak