myprojects-hg-reborn

Fork

#! /usr/bin/env python

import scipy as s
import numpy as n
import pylab as p
import scipy.integrate as si


def simple_coupling(y,t):
    coupling=s.zeros(len(y))

    y_new=y/x_grid
    
    coupling[1:-1]=y_new[2:]-y_new[:-2]
    coupling=-coupling*A/(2.*dx)
    return coupling

A=0.00016358757687

#fix grid in x (careful to avoid zero!)

NN=2000

x_grid=s.linspace(0.01,3.,NN)
print "x_grid is, ", x_grid

dx=x_grid[2]-x_grid[1]

print "dx is, ", dx

#now create a constant initial state

my_sel=s.where((x_grid>=0.4) & (x_grid<=0.5))

# print "x_grid[my_sel] is, ", x_grid[my_sel]

# y0=s.zeros(NN)
# y0[my_sel]=1000.

y0=s.exp(-(x_grid-0.5)**2./0.02)

print "y0 is, ", y0

t=s.linspace(0.,3500., 60)

print "t is, ", t

# y = integrate.odeint(coupling_optim, y0, \
# 	t,printmessg=1,rtol=1e-10,atol=1e-10)


y = si.odeint(simple_coupling, y0, \
	t,printmessg=1, mxstep=600)

#p.save("evolution_of_the_population.dat",y)

print "y[20,:] is, ", y[20,:]

fig = p.figure()
axes = fig.gca()

# p.semilogx(dp,my_initial_state,"ro")
# p.semilogx(dp[selection],grown_state[selection],"kx")
p.plot(x_grid,y0,"ro")
p.plot(x_grid,y[(len(t)-1),:],"b")

p.legend(("Initial state", "10 minutes"))


p.subplots_adjust(bottom=0.2)




p.xlabel(r'$d_p [\mu m]$', fontsize=20)
p.ylabel(r'$n(d_p,t)$', fontsize=20)
p.savefig("evolution_semi_difference.pdf")
 
p.clf()    



print "So far so good"