3. C programming

The following code is needed only if you are using the notebook, because make command treats spaces abd TAB characters differently. We need to be able to reconfigure Jupyter to treat TABS as characters, not replace them with 4 spaces, as it does by default.

In [1]:
%%javascript
IPython.tab_as_tab_everywhere = function(use_tabs) {
    if (use_tabs === undefined) {
        use_tabs = true; 
        }
    // apply setting to all current CodeMirror instances
    IPython.notebook.get_cells().map(
        function(c) {  return c.code_mirror.options.indentWithTabs=use_tabs;  }
        );
    // make sure new CodeMirror instances created in the future also use this setting
    CodeMirror.defaults.indentWithTabs=use_tabs;
    };

IPython.tab_as_tab_everywhere()

Gaussian packet evolution: a C project

More homework! This is the second part of the homework, which is to try to reproduce the functionality of an old Fortran program. We do not know Fortran and little of C syntax, but hopefully the computational content, i.e. the algorithmic core, will be self-explanatory. The task is to "port", or convert this old code to C, and learn some basic C programming in the process.

In [2]:
%%bash

## first, a refresher: what does this bash script do??  Add a few comment lines as appropriate

if [ ! -d ~/5P10 ]; then
  mkdir ~/5P10
fi
cd ~/5P10
rm -rf packet
mkdir packet
In [3]:
%cd ~/5P10/packet
/home/esternin/5P10/packet
In [4]:
%%bash
pwd
ls -l
/home/esternin/5P10/packet
total 0

First, let's review the old Fortran code to see how it works and what its output looks like.

In [5]:
%%file packet.f
C packet.f -- propagation of a Gaussian wave  packet through a barrier in 1D
C	x = coordinate (in units of a, the barrier is from -a to +a)
C	P = wavefunction  (not normalized)
C	E = energy of the wavefunction (in units of V, the barrier height)
C
C	Written by:	E.Sternin
C	Completed:	30.03.93
C=======================================================================	
	real*4		x(601),x0,p0,sp,p_now,p_step,Pi,A0,t
	complex*8	P(601),I,A
	integer		j,Nargs,Nwaves
	character*80 	buffer

	Nwaves = 20	! add up (2*Nwaves+1) eigenwaves to form a packet
	x0 = -30.	! start the packet far away from the barrier
	p0 = 0.5	! keep incident energy low, E=0.5**2 < 1 (units of V)
	sp = 0.1	! p0 +/- 4 sigma_p > 0,	all eigenwaves move to the right
	 
	Nargs = iargc() 	! arguments supplied on the command line ?
	if (Nargs .lt. 1) then 
	  write(*,'(A,$)') '  Error: need time t: '
	  read(*,*) t
	else
	  call getarg(1,buffer)
	  read(buffer,*) t
	end if

	I = cmplx(0.,1.)		!complex i
	Pi = 2.*acos(0.)
	A0 = 1./sqrt( sqrt(2.*Pi) * sp ) /float (2*Nwaves+1)

	do j = 1,601
	  x(j) = -75 + (j-1)*0.25
	  P(j) = cmplx(0.,0.)
	end do

	p_step = 4.*sp/float(Nwaves)	! calculate over +/- 4*sigma_p
	do j = -Nwaves,Nwaves		! in Nwaves steps
	  p_now = p0 + j * p_step 
	  E = p_now**2
	  A = A0*exp( - ( (p_now-p0)/(2*sp) )**2 - I*(p_now*x0 + E*t) )
	  call add_wave(x,P,A,E)
	end do

C..output the saved values as is, without normalization
  10	do j=1,	601
	  write(*,*) x(j),abs(P(j))**2,real(P(j)),aimag(P(j))
	end do
	end

C=======================================================================	
	subroutine add_wave(x,P,A,E)
C..adds the specified wave to the the supplied complex storage array P
C  the values of Psi are only calculated at points given in real array x 
	real*4 		x(601),E
	complex*8	P(601),I,A,k1,K2,C1,b,c,d,f

	I = cmplx(0.,1.)	!complex i

	k1 = sqrt(cmplx(E,0.))
	K2 = sqrt(cmplx(E-1.,0.))

	C1 = -2*k1*K2-k1**2-2*k1*K2*exp(4*I*K2)+exp(4*I*K2)*k1**2
     *	-K2**2+K2**2*exp(4*I*K2)

	b = (k1**2-K2**2)/ C1 *(exp(2*I*(2*K2-k1))-exp(-2*I*k1))
	c = -2*k1*(K2+k1)/ C1 *exp(I*(K2-k1))
	d = (-2*K2+2*k1)*k1/ C1 *exp(I*(-k1+3*K2))
	f = -4*k1*K2/ C1 *exp(2*I*(K2-k1))

	do j = 1,296		! left of the barrier, x < -1
	  P(j) = P(j) + A * ( exp(I*k1*x(j)) + b*exp(-I*k1*x(j)) )
	end do
	do j=297,304		! inside the barrier, -1 < x < 1
	  P(j) = P(j) + A * ( c*exp(I*K2*x(j)) + d*exp(-I*K2*x(j)) )
	end do
	do j=305,601		! to the right of the barrier, x > 1
	  P(j) = P(j) + A * f*exp(I*k1*x(j))
	end do
	  
 	return
	end
Writing packet.f
In [18]:
%%bash
gfortran -o packet packet.f
./packet 40
  -75.0000000       1.96023257E-06   8.77547078E-04   1.09093706E-03
  -74.7500000       1.82045062E-06   1.01503159E-03   8.88910319E-04
  -74.5000000       1.73418891E-06   1.13934174E-03   6.60370453E-04
  -74.2500000       1.72351565E-06   1.24439807E-03   4.18317009E-04
  -74.0000000       1.78586924E-06   1.32473465E-03   1.75918191E-04
  -73.7500000       1.89733987E-06   1.37641025E-03  -5.32418489E-05
  -73.5000000       2.01731950E-06   1.39714230E-03  -2.55563966E-04
  -73.2500000       2.09439304E-06   1.38478726E-03  -4.20425145E-04
  -73.0000000       2.08815527E-06   1.34157971E-03  -5.36953798E-04
  -72.7500000       1.96922588E-06   1.26916065E-03  -5.98713057E-04
  -72.5000000       1.73377737E-06   1.17202860E-03  -6.00105326E-04
  -72.2500000       1.40628003E-06   1.05583435E-03  -5.39901666E-04
  -72.0000000       1.03724994E-06   9.27589601E-04  -4.20508557E-04
  -71.7500000       6.92532581E-07   7.95061758E-04  -2.45783187E-04
  -71.5000000       4.43797006E-07   6.65725674E-04  -2.46243435E-05
  -71.2500000       3.53931824E-07   5.47598756E-04   2.32524064E-04
  -71.0000000       4.64492956E-07   4.48164472E-04   5.13460371E-04
  -70.7500000       7.85790803E-07   3.73223447E-04   8.04049196E-04
  -70.5000000       1.29298166E-06   3.28127062E-04   1.08872144E-03
  -70.2500000       1.93195683E-06   3.15745245E-04   1.35361066E-03
  -70.0000000       2.62589083E-06   3.38269514E-04   1.58476015E-03
  -69.7500000       3.28651595E-06   3.95180687E-04   1.76927901E-03
  -69.5000000       3.83478164E-06   4.83024662E-04   1.89775368E-03
  -69.2500000       4.21086088E-06   5.98888146E-04   1.96270063E-03
  -69.0000000       4.38410598E-06   7.35657639E-04   1.96033507E-03
  -68.7500000       4.35984521E-06   8.86410533E-04   1.89053488E-03
  -68.5000000       4.17116462E-06   1.04239804E-03   1.75629475E-03
  -68.2500000       3.87792215E-06   1.19528465E-03   1.56499748E-03
  -68.0000000       3.54403505E-06   1.33645989E-03   1.32586202E-03
  -67.7500000       3.23056088E-06   1.45724276E-03   1.05214270E-03
  -67.5000000       2.97883844E-06   1.55024452E-03   7.58670038E-04
  -67.2500000       2.80352128E-06   1.60967337E-03   4.60947951E-04
  -67.0000000       2.69476050E-06   1.63207762E-03   1.76304573E-04
  -66.7500000       2.61299533E-06   1.61450752E-03  -7.97556131E-05
  -66.5000000       2.51203733E-06   1.55792257E-03  -2.91400909E-04
  -66.2500000       2.34614959E-06   1.46510801E-03  -4.46775230E-04
  -66.0000000       2.08429469E-06   1.34065363E-03  -5.35670319E-04
  -65.7500000       1.72119576E-06   1.19045167E-03  -5.51380508E-04
  -65.5000000       1.28981821E-06   1.02388021E-03  -4.91413812E-04
  -65.2500000       8.50275796E-07   8.49914155E-04  -3.57661454E-04
  -65.0000000       4.84994189E-07   6.79052668E-04  -1.54537032E-04
  -64.7500000       2.82709209E-07   5.20876376E-04   1.06756575E-04
  -64.5000000       3.20579090E-07   3.85980820E-04   4.14243783E-04
  -64.2500000       6.45942976E-07   2.82163615E-04   7.52546766E-04
  -64.0000000       1.26368809E-06   2.17078952E-04   1.10297999E-03
  -63.7500000       2.13198837E-06   1.95043162E-04   1.44704757E-03
  -63.5000000       3.16631758E-06   2.17700610E-04   1.76604756E-03
  -63.2500000       4.24724749E-06   2.85472372E-04   2.04101764E-03
  -63.0000000       5.24817187E-06   3.94371396E-04   2.25668866E-03
  -62.7500000       6.04795923E-06   5.39795437E-04   2.39928742E-03
  -62.5000000       6.55504118E-06   7.12396111E-04   2.45917309E-03
  -62.2500000       6.72689112E-06   9.02541215E-04   2.43152422E-03
  -62.0000000       6.56544489E-06   1.09797518E-03   2.31514475E-03
  -61.7500000       6.12861822E-06   1.28629222E-03   2.11520004E-03
  -61.5000000       5.50341974E-06   1.45479967E-03   1.84037443E-03
  -61.2500000       4.79634036E-06   1.59180793E-03   1.50415674E-03
  -61.0000000       4.10773055E-06   1.68611656E-03   1.12460728E-03
  -60.7500000       3.51199446E-06   1.72967336E-03   7.21265911E-04
  -60.5000000       3.04656032E-06   1.71629013E-03   3.17660917E-04
  -60.2500000       2.70658597E-06   1.64393231E-03  -6.38172496E-05
  -60.0000000       2.44762691E-06   1.51251513E-03  -3.99906014E-04
  -59.7500000       2.21115170E-06   1.32758426E-03  -6.69829722E-04
  -59.5000000       1.93523533E-06   1.09681441E-03  -8.55706516E-04
  -59.2500000       1.58132946E-06   8.31471640E-04  -9.43389896E-04
  -59.0000000       1.15117166E-06   5.46015042E-04  -9.23601212E-04
  -58.7500000       6.96210464E-07   2.55807769E-04  -7.94212101E-04
  -58.5000000       3.10616684E-07  -2.11259467E-05  -5.56929444E-04
  -58.2500000       1.20390737E-07  -2.67395400E-04  -2.21111812E-04
  -58.0000000       2.56541142E-07  -4.66421829E-04   1.97463436E-04
  -57.7500000       8.25106667E-07  -6.03376422E-04   6.79001911E-04
  -57.5000000       1.87918408E-06  -6.65937783E-04   1.19821157E-03
  -57.2500000       3.39837993E-06  -6.45499444E-04   1.72676297E-03
  -57.0000000       5.28219744E-06  -5.38554392E-04   2.23431340E-03
  -56.7500000       7.35717913E-06  -3.46145709E-04   2.69023469E-03
  -56.5000000       9.40185237E-06  -7.51864864E-05   3.06532206E-03
  -56.2500000       1.11809504E-05   2.63452530E-04   3.33339814E-03
  -56.0000000       1.24820490E-05   6.52841758E-04   3.47215310E-03
  -55.7500000       1.31673269E-05   1.07204251E-03   3.46670626E-03
  -55.5000000       1.31818224E-05   1.49678183E-03   3.30778863E-03
  -55.2500000       1.25843872E-05   1.90045475E-03   2.99543981E-03
  -55.0000000       1.15273106E-05   2.25543324E-03   2.53778091E-03
  -54.7500000       1.02349368E-05   2.53525679E-03   1.95125840E-03
  -54.5000000       8.96054371E-06   2.71489052E-03   1.26091787E-03
  -54.2500000       7.94261359E-06   2.77374568E-03   4.98947629E-04
  -54.0000000       7.35255207E-06   2.69523030E-03  -2.97128689E-04
  -53.7500000       7.28275381E-06   2.47116061E-03  -1.08449056E-03
  -53.5000000       7.71300893E-06   2.09917175E-03  -1.81837473E-03
  -53.2500000       8.53753500E-06   1.58579578E-03  -2.45413673E-03
  -53.0000000       9.59492263E-06   9.45997075E-04  -2.94957822E-03
  -52.7500000       1.07135056E-05   2.03106902E-04  -3.26684164E-03
  -52.5000000       1.17637628E-05  -6.11433934E-04  -3.37489438E-03
  -52.2500000       1.27123058E-05  -1.46040926E-03  -3.25261592E-03
  -52.0000000       1.36412191E-05  -2.30117701E-03  -2.88891047E-03
  -51.7500000       1.47498858E-05  -3.08695436E-03  -2.28486280E-03
  -51.5000000       1.63353143E-05  -3.77088925E-03  -1.45454728E-03
  -51.2500000       1.87299502E-05  -4.30688309E-03  -4.25099803E-04
  -51.0000000       2.22337785E-05  -4.65292996E-03   7.64213735E-04
  -50.7500000       2.70391229E-05  -4.77324473E-03   2.06282735E-03
  -50.5000000       3.31695010E-05  -4.64084372E-03   3.41058243E-03
  -50.2500000       4.04509083E-05  -4.23978781E-03   4.74079186E-03
  -50.0000000       4.84937918E-05  -3.56642436E-03   5.98117150E-03
  -49.7500000       5.67711140E-05  -2.63119512E-03   7.06030661E-03
  -49.5000000       6.46740809E-05  -1.45877257E-03   7.90860690E-03
  -49.2500000       7.16222858E-05  -8.82214517E-05   8.46253522E-03
  -49.0000000       7.71898267E-05   1.42793136E-03   8.66895821E-03
  -48.7500000       8.11878417E-05   3.02483817E-03   8.48753192E-03
  -48.5000000       8.37276530E-05   4.62769996E-03   7.89379794E-03
  -48.2500000       8.52661542E-05   6.15607155E-03   6.88250922E-03
  -48.0000000       8.65038237E-05   7.52407406E-03   5.46736922E-03
  -47.7500000       8.83546600E-05   8.64779390E-03   3.68379219E-03
  -47.5000000       9.17651269E-05   9.44703538E-03   1.58702326E-03
  -47.2500000       9.75877483E-05   9.85029712E-03  -7.47931248E-04
  -47.0000000       1.06456144E-04   9.79940034E-03  -3.22922552E-03
  -46.7500000       1.18670869E-04   9.25225951E-03  -5.75035345E-03
  -46.5000000       1.34186004E-04   8.18734709E-03  -8.19471478E-03
  -46.2500000       1.52631736E-04   6.60495553E-03  -1.04406076E-02
  -46.0000000       1.73423046E-04   4.53073671E-03  -1.23650907E-02
  -45.7500000       1.95925590E-04   2.01440044E-03  -1.38516352E-02
  -45.5000000       2.19618625E-04  -8.67994619E-04  -1.47940936E-02
  -45.2500000       2.44244642E-04  -4.01895726E-03  -1.51027357E-02
  -45.0000000       2.69948883E-04  -7.31897354E-03  -1.47099122E-02
  -44.7500000       2.97317951E-04  -1.06324060E-02  -1.35746058E-02
  -44.5000000       3.27306130E-04  -1.38114765E-02  -1.16854291E-02
  -44.2500000       3.61102109E-04  -1.67018082E-02  -9.06375702E-03
  -44.0000000       3.99954821E-04  -1.91496704E-02  -5.76584507E-03
  -43.7500000       4.44861274E-04  -2.10077558E-02  -1.88028580E-03
  -43.5000000       4.96436493E-04  -2.21435502E-02   2.46975827E-03
  -43.2500000       5.54633385E-04  -2.24444978E-02   7.13287760E-03
  -43.0000000       6.18779508E-04  -2.18262710E-02   1.19328713E-02
  -42.7500000       6.87542197E-04  -2.02368461E-02   1.66736990E-02
  -42.5000000       7.59166374E-04  -1.76626258E-02   2.11470556E-02
  -42.2500000       8.31750222E-04  -1.41301081E-02   2.51414049E-02
  -42.0000000       9.03629640E-04  -9.71046370E-03   2.84488406E-02
  -41.7500000       9.73681279E-04  -4.51691356E-03   3.08752116E-02
  -41.5000000       1.04165531E-03   1.29624316E-03   3.22486460E-02
  -41.2500000       1.10842264E-03   7.53587764E-03   3.24288942E-02
  -41.0000000       1.17590069E-03   1.39767462E-02   3.13137546E-02
  -40.7500000       1.24702416E-03   2.03675032E-02   2.88476888E-02
  -40.5000000       1.32536935E-03   2.64405254E-02   2.50253454E-02
  -40.2500000       1.41476106E-03   3.19203436E-02   1.98960491E-02
  -40.0000000       1.51880877E-03   3.65348384E-02   1.35651864E-02
  -39.7500000       1.64052797E-03   4.00270075E-02   6.19407417E-03
  -39.5000000       1.78196223E-03   4.21656929E-02  -2.00418197E-03
  -39.2500000       1.94406696E-03   4.27564979E-02  -1.07679572E-02
  -39.0000000       2.12683622E-03   4.16527316E-02  -1.97961163E-02
  -38.7500000       2.32949387E-03   3.87631059E-02  -2.87561417E-02
  -38.5000000       2.55105319E-03   3.40603143E-02  -3.72954160E-02
  -38.2500000       2.79075885E-03   2.75851320E-02  -4.50535156E-02
  -38.0000000       3.04869516E-03   1.94496233E-02  -5.16759828E-02
  -37.7500000       3.32612381E-03   9.83738527E-03  -5.68273664E-02
  -37.5000000       3.62580363E-03  -9.99222975E-04  -6.02063537E-02
  -37.2500000       3.95175675E-03  -1.27440859E-02  -6.15576580E-02
  -37.0000000       4.30907169E-03  -2.50263419E-02  -6.06856942E-02
  -36.7500000       4.70315013E-03  -3.74297015E-02  -5.74644916E-02
  -36.5000000       5.13906637E-03  -4.95071150E-02  -5.18469997E-02
  -36.2500000       5.62049542E-03  -6.07937984E-02  -4.38703634E-02
  -36.0000000       6.14928734E-03  -7.08252192E-02  -3.36611867E-02
  -35.7500000       6.72455505E-03  -7.91526288E-02  -2.14340165E-02
  -35.5000000       7.34266499E-03  -8.53614360E-02  -7.48931477E-03
  -35.2500000       7.99756218E-03  -8.90889466E-02   7.79240625E-03
  -35.0000000       8.68107937E-03  -9.00388882E-02   2.39599366E-02
  -34.7500000       9.38394945E-03  -8.79963562E-02   4.05042097E-02
  -34.5000000       1.00970156E-02  -8.28404799E-02   5.68724088E-02
  -34.2500000       1.08124306E-02  -7.45537356E-02   7.24856555E-02
  -34.0000000       1.15246875E-02  -6.32275492E-02   8.67580846E-02
  -33.7500000       1.22313322E-02  -4.90646362E-02   9.91160572E-02
  -33.5000000       1.29336873E-02  -3.23789082E-02  0.109019697    
  -33.2500000       1.36364633E-02  -1.35887824E-02  0.115981929    
  -33.0000000       1.43475989E-02   6.79147011E-03  0.119588777    
  -32.7500000       1.50774391E-02   2.81644873E-02  0.119516529    
  -32.5000000       1.58378147E-02   4.98665795E-02  0.115547128    
  -32.2500000       1.66410562E-02   7.11870566E-02  0.107580006    
  -32.0000000       1.74995102E-02   9.13899839E-02   9.56419334E-02
  -31.7500000       1.84248574E-02  0.109737672       7.98905417E-02
  -31.5000000       1.94284152E-02  0.125515372       6.06160499E-02
  -31.2500000       2.05216613E-02  0.138056397       3.82373072E-02
  -31.0000000       2.17168089E-02  0.146765426       1.32935783E-02
  -30.7500000       2.30280850E-02  0.151142314      -1.35678239E-02
  -30.5000000       2.44722031E-02  0.150801063      -4.16081809E-02
  -30.2500000       2.60691196E-02  0.145488039      -7.00167865E-02
  -30.0000000       2.78420150E-02  0.135095716      -9.79344919E-02
  -29.7500000       2.98161153E-02  0.119670928     -0.124478847    
  -29.5000000       3.20172086E-02   9.94206443E-02 -0.148770779    
  -29.2500000       3.44693512E-02   7.47105330E-02 -0.169963777    
  -29.0000000       3.71915400E-02   4.60601375E-02 -0.187269866    
  -28.7500000       4.01945785E-02   1.41317565E-02 -0.199987173    
  -28.5000000       4.34780344E-02  -2.02847160E-02 -0.207524851    
  -28.2500000       4.70276065E-02  -5.62914461E-02 -0.209425122    
  -28.0000000       5.08133098E-02  -9.29070786E-02 -0.205381557    
  -27.7500000       5.47894165E-02 -0.129093677     -0.195254281    
  -27.5000000       5.88950552E-02 -0.163785622     -0.179079115    
  -27.2500000       6.30571246E-02 -0.195921570     -0.157072783    
  -27.0000000       6.71939254E-02 -0.224476039     -0.129631907    
  -26.7500000       7.12202266E-02 -0.248490781      -9.73270535E-02
  -26.5000000       7.50526562E-02 -0.267104775      -6.08909056E-02
  -26.2500000       7.86157101E-02 -0.279582143      -2.12023389E-02
  -26.0000000       8.18473473E-02 -0.285337299       2.07357500E-02
  -25.7500000       8.47036093E-02 -0.283954978       6.38214201E-02
  -25.5000000       8.71627629E-02 -0.275206804      0.106883004    
  -25.2500000       8.92279223E-02 -0.259061873      0.148710683    
  -25.0000000       9.09286290E-02 -0.235692784      0.188089162    
  -24.7500000       9.23208743E-02 -0.205473825      0.223833382    
  -24.5000000       9.34863836E-02 -0.168976292      0.254820317    
  -24.2500000       9.45309997E-02 -0.126955181      0.280023903    
  -24.0000000       9.55815017E-02  -8.03322792E-02  0.298543513    
  -23.7500000       9.67830643E-02  -3.01735159E-02  0.309633046    
  -23.5000000       9.82951000E-02   2.23375410E-02  0.312723726    
  -23.2500000      0.100286402       7.59295449E-02  0.307442844    
  -23.0000000      0.102930218      0.129277900      0.293628067    
  -22.7500000      0.106397383      0.181037560      0.271335185    
  -22.5000000      0.110849135      0.229880571      0.240840301    
  -22.2500000      0.116428502      0.274530202      0.202636763    
  -22.0000000      0.123250365      0.313795745      0.157424912    
  -21.7500000      0.131391481      0.346604764      0.106097236    
  -21.5000000      0.140879929      0.372032285       4.97182272E-02
  -21.2500000      0.151684910      0.389325887      -1.05014276E-02
  -21.0000000      0.163708612      0.397926509      -7.32332170E-02
  -20.7500000      0.176780686      0.397485137     -0.137062952    
  -20.5000000      0.190655395      0.387872219     -0.200525656    
  -20.2500000      0.205013812      0.369182646     -0.262141168    
  -20.0000000      0.219471708      0.341735393     -0.320450664    
  -19.7500000      0.233590782      0.306065470     -0.374051750    
  -19.5000000      0.246898264      0.262913287     -0.421633571    
  -19.2500000      0.258908153      0.213205934     -0.462007999    
  -19.0000000      0.269148082      0.158035666     -0.494138479    
  -18.7500000      0.277188301       9.86343697E-02 -0.517164946    
  -18.5000000      0.282670110       3.63434292E-02 -0.530423641    
  -18.2500000      0.285334557      -2.74175666E-02 -0.533463061    
  -18.0000000      0.285046995      -9.11845043E-02 -0.526053607    
  -17.7500000      0.281816214     -0.153482959     -0.508192062    
  -17.5000000      0.275807947     -0.212862998     -0.480101347    
  -17.2500000      0.267349601     -0.267932624     -0.442223638    
  -17.0000000      0.256926596     -0.317389399     -0.395209521    
  -16.7500000      0.245169237     -0.360050529     -0.339901239    
  -16.5000000      0.232831746     -0.394878566     -0.277313322    
  -16.2500000      0.220762447     -0.421005070     -0.208607733    
  -16.0000000      0.209867477     -0.437749028     -0.135067612    
  -15.7500000      0.201068863     -0.444631368      -5.80672622E-02
  -15.5000000      0.195259273     -0.441384226       2.09580418E-02
  -15.2500000      0.193256199     -0.427956969      0.100543648    
  -15.0000000      0.195754766     -0.404514939      0.179227337    
  -14.7500000      0.203285381     -0.371435642      0.255579591    
  -14.5000000      0.216175333     -0.329299301      0.328233600    
  -14.2500000      0.234517634     -0.278875738      0.395911574    
  -14.0000000      0.258148938     -0.221106887      0.457450181    
  -13.7500000      0.286637366     -0.157086968      0.511821330    
  -13.5000000      0.319282383      -8.80380198E-02  0.558150232    
  -13.2500000      0.355127424      -1.52855385E-02  0.595729589    
  -13.0000000      0.392985851       5.97697571E-02  0.624029994    
  -12.7500000      0.431478918      0.135679722      0.642705142    
  -12.5000000      0.469087094      0.210978046      0.651594460    
  -12.2500000      0.504212439      0.284208447      0.650721133    
  -12.0000000      0.535248458      0.353953332      0.640285492    
  -11.7500000      0.560658395      0.418859571      0.620656967    
  -11.5000000      0.579053819      0.477663875      0.592360556    
  -11.2500000      0.589274645      0.529216051      0.556062043    
  -11.0000000      0.590462565      0.572498322      0.512550712    
  -10.7500000      0.582124531      0.606643021      0.462718904    
  -10.5000000      0.564184606      0.630946636      0.407542557    
  -10.2500000      0.537014902      0.644880176      0.348058075    
  -10.0000000      0.501449108      0.648096800      0.285341263    
  -9.75000000      0.458772004      0.640436172      0.220484659    
  -9.50000000      0.410683244      0.621923983      0.154576093    
  -9.25000000      0.359239548      0.592769504       8.86783823E-02
  -9.00000000      0.306772977      0.553358912       2.38095820E-02
  -8.75000000      0.255792588      0.504247785      -3.90739515E-02
  -8.50000000      0.208866581      0.446146846      -9.90937650E-02
  -8.25000000      0.168500155      0.379909933     -0.155462518    
  -8.00000000      0.137006804      0.306517273     -0.207494527    
  -7.75000000      0.116383560      0.227056608     -0.254615128    
  -7.50000000      0.108197428      0.142704770     -0.296365947    
  -7.25000000      0.113488071       5.47072664E-02 -0.332408160    
  -7.00000000      0.132693723      -3.56425345E-02 -0.362523556    
  -6.75000000      0.165604278     -0.127024367     -0.386612326    
  -6.50000000      0.211346462     -0.218111813     -0.404689640    
  -6.25000000      0.268400699     -0.307592630     -0.416878223    
  -6.00000000      0.334652424     -0.394188046     -0.423400760    
  -5.75000000      0.407474846     -0.476671308     -0.424569577    
  -5.50000000      0.483839244     -0.553883970     -0.420775205    
  -5.25000000      0.560449958     -0.624752462     -0.412473440    
  -5.00000000      0.633894145     -0.688300073     -0.400171399    
  -4.75000000      0.700804830     -0.743660212     -0.384414226    
  -4.50000000      0.758022249     -0.790085316     -0.365769655    
  -4.25000000      0.802751303     -0.826954842     -0.344814390    
  -4.00000000      0.832703292     -0.853781283     -0.322119206    
  -3.75000000      0.846216321     -0.870213330     -0.298236668    
  -3.50000000      0.842348874     -0.876038551     -0.273688346    
  -3.25000000      0.820938110     -0.871183276     -0.248953462    
  -3.00000000      0.782622397     -0.855710447     -0.224459350    
  -2.75000000      0.728825748     -0.829816759     -0.200573817    
  -2.50000000      0.661704242     -0.793828130     -0.177598387    
  -2.25000000      0.584054708     -0.748192668     -0.155763954    
  -2.00000000      0.499194533     -0.693475068     -0.135228887    
  -1.75000000      0.410810202     -0.630346000     -0.116078191    
  -1.50000000      0.322792500     -0.559575438      -9.83251929E-02
  -1.25000000      0.239054367     -0.482021004      -8.19153488E-02
  -1.00000000      0.163349345     -0.398617983      -6.67311698E-02
 -0.750000000      0.106638193     -0.321966916      -5.45481965E-02
 -0.500000000       7.05969557E-02 -0.261571050      -4.66640629E-02
 -0.250000000       4.77470905E-02 -0.214356616      -4.24067602E-02
   0.00000000       3.33616547E-02 -0.177911237      -4.13430445E-02
  0.250000000       2.44787242E-02 -0.150361016      -4.32468429E-02
  0.500000000       1.92833506E-02 -0.130275324      -4.80800532E-02
  0.750000000       1.67281050E-02 -0.116592705      -5.59843406E-02
   1.00000000       1.63136702E-02 -0.108565912      -6.72838166E-02
   1.25000000       1.69207063E-02 -0.102376230      -8.02484378E-02
   1.50000000       1.75367892E-02  -9.45702568E-02  -9.26998109E-02
   1.75000000       1.81611106E-02  -8.51789638E-02 -0.104430154    
   2.00000000       1.87927652E-02  -7.42605627E-02 -0.115230784    
   2.25000000       1.94308497E-02  -6.19015098E-02 -0.124896169    
   2.50000000       2.00744048E-02  -4.82166074E-02 -0.133227497    
   2.75000000       2.07224004E-02  -3.33488137E-02 -0.140036628    
   3.00000000       2.13738922E-02  -1.74683388E-02 -0.145150781    
   3.25000000       2.20277868E-02  -7.71210820E-04 -0.148415610    
   3.50000000       2.26830319E-02   1.65228322E-02 -0.149699792    
   3.75000000       2.33385805E-02   3.41729783E-02 -0.148898587    
   4.00000000       2.39933450E-02   5.19205667E-02 -0.145936966    
   4.25000000       2.46462654E-02   6.94929063E-02 -0.140772864    
   4.50000000       2.52962392E-02   8.66073593E-02 -0.133399412    
   4.75000000       2.59422213E-02  0.102976345     -0.123847060    
   5.00000000       2.65831538E-02  0.118311986     -0.112184785    
   5.25000000       2.72179805E-02  0.132331669      -9.85206142E-02
   5.50000000       2.78457049E-02  0.144763365      -8.30016583E-02
   5.75000000       2.84652840E-02  0.155350968      -6.58130646E-02
   6.00000000       2.90757511E-02  0.163860142      -4.71763238E-02
   6.25000000       2.96761319E-02  0.170083150      -2.73469239E-02
   6.50000000       3.02655026E-02  0.173844174      -6.61101239E-03
   6.75000000       3.08429208E-02  0.175003663       1.47186136E-02
   7.00000000       3.14075276E-02  0.173462734       3.63071337E-02
   7.25000000       3.19584385E-02  0.169166327       5.78030497E-02
   7.50000000       3.24948132E-02  0.162106231       7.88440332E-02
   7.75000000       3.30158845E-02  0.152323157       9.90633070E-02
   8.00000000       3.35208140E-02  0.139907509      0.118096165    
   8.25000000       3.40088643E-02  0.124999985      0.135587126    
   8.50000000       3.44793163E-02  0.107790716      0.151196808    
   8.75000000       3.49314287E-02   8.85174349E-02  0.164608911    
   9.00000000       3.53645310E-02   6.74628466E-02  0.175537154    
   9.25000000       3.57779264E-02   4.49507721E-02  0.183731750    
   9.50000000       3.61710414E-02   2.13416759E-02  0.188985646    
   9.75000000       3.65432091E-02  -2.97324592E-03  0.191139653    
   10.0000000       3.68938185E-02  -2.75782626E-02  0.190087497    
   10.2500000       3.72223109E-02  -5.20405248E-02  0.185779691    
   10.5000000       3.75281312E-02  -7.59176388E-02  0.178226382    
   10.7500000       3.78107764E-02  -9.87662971E-02  0.167499229    
   11.0000000       3.80696766E-02 -0.120150708      0.153731868    
   11.2500000       3.83044183E-02 -0.139651760      0.137119666    
   11.5000000       3.85145098E-02 -0.156875610      0.117917560    
   11.7500000       3.86995003E-02 -0.171462551       9.64369774E-02
   12.0000000       3.88590544E-02 -0.183095634       7.30413720E-02
   12.2500000       3.89927216E-02 -0.191507757       4.81404029E-02
   12.5000000       3.91002409E-02 -0.196489573       2.21830513E-02
   12.7500000       3.91812734E-02 -0.197894782      -4.35046153E-03
   13.0000000       3.92355770E-02 -0.195645571      -3.09577901E-02
   13.2500000       3.92628834E-02 -0.189735904      -5.71241416E-02
   13.5000000       3.92630734E-02 -0.180234164      -8.23329613E-02
   13.7500000       3.92359607E-02 -0.167283267     -0.106076717    
   14.0000000       3.91815081E-02 -0.151100218     -0.127868026    
   14.2500000       3.90996300E-02 -0.131972954     -0.147250712    
   14.5000000       3.89903784E-02 -0.110256314     -0.163810641    
   14.7500000       3.88537943E-02  -8.63659829E-02 -0.177185535    
   15.0000000       3.86900119E-02  -6.07706420E-02 -0.187074691    
   15.2500000       3.84992287E-02  -3.39830369E-02 -0.193246946    
   15.5000000       3.82816680E-02  -6.54960843E-03 -0.195547357    
   15.7500000       3.80376503E-02   2.09612399E-02 -0.193902746    
   16.0000000       3.77674885E-02   4.79713529E-02 -0.188324824    
   16.2500000       3.74716371E-02   7.39059523E-02 -0.178912118    
   16.5000000       3.71505991E-02   9.82069075E-02 -0.165849343    
   16.7500000       3.68048996E-02  0.120346002     -0.149404615    
   17.0000000       3.64351235E-02  0.139837801     -0.129925027    
   17.2500000       3.60419489E-02  0.156251878     -0.107829943    
   17.5000000       3.56260724E-02  0.169224024      -8.36020410E-02
   17.7500000       3.51882987E-02  0.178465948      -5.77771701E-02
   18.0000000       3.47294174E-02  0.183773309      -3.09320241E-02
   18.2500000       3.42502892E-02  0.185031936      -3.67091759E-03
   18.5000000       3.37518603E-02  0.182221964       2.33883467E-02
   18.7500000       3.32350954E-02  0.175419420       4.96298186E-02
   19.0000000       3.27009447E-02  0.164795294       7.44543746E-02
   19.2500000       3.21504772E-02  0.150612742       9.72948074E-02
   19.5000000       3.15847397E-02  0.133221060      0.117630303    
   19.7500000       3.10048386E-02  0.113047987      0.134999976    
   20.0000000       3.04118637E-02   9.05891880E-02  0.149014980    
   20.2500000       2.98069343E-02   6.63963482E-02  0.159368947    
   20.5000000       2.91912276E-02   4.10635658E-02  0.165846348    
   20.7500000       2.85658613E-02   1.52118383E-02  0.168328434    
   21.0000000       2.79319920E-02  -1.05265593E-02  0.166796833    
   21.2500000       2.72907782E-02  -3.55247296E-02  0.161334351    
   21.5000000       2.66433768E-02  -5.91774285E-02  0.152123004    
   21.7500000       2.59909071E-02  -8.09171647E-02  0.139439315    
   22.0000000       2.53345203E-02 -0.100229234      0.123647168    
   22.2500000       2.46753097E-02 -0.116665415      0.105187878    
   22.5000000       2.40143724E-02 -0.129855871       8.45684558E-02
   22.7500000       2.33527794E-02 -0.139519006       6.23476058E-02
   23.0000000       2.26915646E-02 -0.145468727       3.91205288E-02
   23.2500000       2.20317636E-02 -0.147619247       1.55023299E-02
   23.5000000       2.13743541E-02 -0.145986721      -7.88856391E-03
   23.7500000       2.07202565E-02 -0.140688226      -3.04479599E-02
   24.0000000       2.00704020E-02 -0.131937996      -5.16020283E-02
   24.2500000       1.94257032E-02 -0.120040581      -7.08234385E-02
   24.5000000       1.87869519E-02 -0.105381221      -8.76455903E-02
   24.7500000       1.81549713E-02  -8.84147659E-02 -0.101674967    
   25.0000000       1.75304953E-02  -6.96514621E-02 -0.112601824    
   25.2500000       1.69142596E-02  -4.96424399E-02 -0.120207690    
   25.5000000       1.63069330E-02  -2.89627798E-02 -0.124370784    
   25.7500000       1.57091450E-02  -8.19526892E-03 -0.125067905    
   26.0000000       1.51214655E-02   1.20870015E-02 -0.122373894    
   26.2500000       1.45444535E-02   3.13367359E-02 -0.116457984    
   26.5000000       1.39785921E-02   4.90483828E-02 -0.107577175    
   26.7500000       1.34243481E-02   6.47718385E-02  -9.60674584E-02
   27.0000000       1.28821200E-02   7.81249329E-02  -8.23323429E-02
   27.2500000       1.23522701E-02   8.88035521E-02  -6.68296218E-02
   27.5000000       1.18351402E-02   9.65890288E-02  -5.00569604E-02
   27.7500000       1.13309966E-02  0.101352930      -3.25358212E-02
   28.0000000       1.08400732E-02  0.103059024      -1.47956591E-02
   28.2500000       1.03625869E-02  0.101762496       2.64230440E-03
   28.5000000       9.89867002E-03   9.76059064E-02   1.92810073E-02
   28.7500000       9.44846403E-03   9.08131301E-02   3.46617773E-02
   29.0000000       9.01204161E-03   8.16802457E-02   4.83774692E-02
   29.2500000       8.58943630E-03   7.05646351E-02   6.00838438E-02
   29.5000000       8.18066113E-03   5.78725114E-02   6.95085153E-02
   29.7500000       7.78569886E-03   4.40449156E-02   7.64574632E-02
   30.0000000       7.40448153E-03   2.95427814E-02   8.08189660E-02
   30.2500000       7.03691412E-03   1.48322899E-02   8.25646296E-02
   30.5000000       6.68288302E-03   3.69874993E-04   8.17480683E-02
   30.7500000       6.34223549E-03  -1.34120928E-02   7.85006434E-02
   31.0000000       6.01479318E-03  -2.61180829E-02   7.30249211E-02
   31.2500000       5.70036517E-03  -3.74010690E-02   6.55860156E-02
   31.5000000       5.39871864E-03  -4.69719134E-02   5.65009527E-02
   31.7500000       5.10962540E-03  -5.46070412E-02   4.61269617E-02
   32.0000000       4.83280607E-03  -6.01531900E-02   3.48482355E-02
   32.2500000       4.56799101E-03  -6.35302737E-02   2.30628457E-02
   32.5000000       4.31489572E-03  -6.47313446E-02   1.11690946E-02
   32.7500000       4.07320121E-03  -6.38200641E-02  -4.47711791E-04
   33.0000000       3.84260272E-03  -6.09263442E-02  -1.14272917E-02
   33.2500000       3.62276752E-03  -5.62392846E-02  -2.14455165E-02
   33.5000000       3.41337430E-03  -4.99988981E-02  -3.02239023E-02
   33.7500000       3.21407593E-03  -4.24858741E-02  -3.75369973E-02
   34.0000000       3.02455202E-03  -3.40108722E-02  -4.32181992E-02
   34.2500000       2.84444378E-03  -2.49024481E-02  -4.71626110E-02
   34.5000000       2.67341733E-03  -1.54955760E-02  -4.93285358E-02
   34.7500000       2.51114182E-03  -6.11969689E-03  -4.97362167E-02
   35.0000000       2.35727243E-03   2.91231833E-03  -4.84643243E-02
   35.2500000       2.21148250E-03   1.13139851E-02  -4.56451140E-02
   35.5000000       2.07344349E-03   1.88343842E-02  -4.14573215E-02
   35.7500000       1.94283889E-03   2.52649616E-02  -3.61181460E-02
   36.0000000       1.81934843E-03   3.04451883E-02  -2.98737157E-02
   36.2500000       1.70266768E-03   3.42659019E-02  -2.29894668E-02
   36.5000000       1.59250724E-03   3.66710722E-02  -1.57397483E-02
   36.7500000       1.48856069E-03   3.76568809E-02  -8.39762110E-03
   37.0000000       1.39055762E-03   3.72700393E-02  -1.22548896E-03
   37.2500000       1.29822362E-03   3.56033407E-02   5.53405704E-03
   37.5000000       1.21128664E-03   3.27901691E-02   1.16658220E-02
   37.7500000       1.12949871E-03   2.89980359E-02   1.69886015E-02
   38.0000000       1.05260930E-03   2.44205855E-02   2.13598758E-02
   38.2500000       9.80380690E-04   1.92693155E-02   2.46794261E-02
   38.5000000       9.12581163E-04   1.37650315E-02   2.68906131E-02
   38.7500000       8.48987489E-04   8.12926143E-03   2.79803965E-02
   39.0000000       7.89391284E-04   2.57591950E-03   2.79777758E-02
   39.2500000       7.33579276E-04  -2.69599282E-03   2.69501563E-02
   39.5000000       6.81357633E-04  -7.50931306E-03   2.49993559E-02
   39.7500000       6.32527866E-04  -1.17140356E-02   2.22555455E-02
   40.0000000       5.86911920E-04  -1.51915830E-02   1.88713465E-02
   40.2500000       5.44333132E-04  -1.78577229E-02   1.50144883E-02
   40.5000000       5.04613039E-04  -1.96636375E-02   1.08606806E-02
   40.7500000       4.67591046E-04  -2.05962434E-02   6.58679148E-03
   41.0000000       4.33109759E-04  -2.06766278E-02   2.36364175E-03
   41.2500000       4.01019206E-04  -1.99573599E-02  -1.65013608E-03
   41.5000000       3.71168629E-04  -1.85185429E-02  -5.31339552E-03
   41.7500000       3.43424355E-04  -1.64636020E-02  -8.50729924E-03
   42.0000000       3.17649014E-04  -1.39135094E-02  -1.11383684E-02
   42.2500000       2.93717079E-04  -1.10013243E-02  -1.31410789E-02
   42.5000000       2.71510566E-04  -7.86614418E-03  -1.44787552E-02
   42.7500000       2.50909536E-04  -4.64736158E-03  -1.51430368E-02
   43.0000000       2.31807731E-04  -1.47880171E-03  -1.51532469E-02
   43.2500000       2.14101543E-04   1.51587930E-03  -1.45534761E-02
   43.5000000       1.97694157E-04   4.22786828E-03  -1.34096714E-02
   43.7500000       1.82495016E-04   6.56684255E-03  -1.18055744E-02
   44.0000000       1.68417479E-04   8.46317131E-03  -9.83830355E-03
   44.2500000       1.55378846E-04   9.86974686E-03  -7.61360303E-03
   44.5000000       1.43307712E-04   1.07628535E-02  -5.24105877E-03
   44.7500000       1.32132191E-04   1.11412527E-02  -2.82925414E-03
   45.0000000       1.21787256E-04   1.10252183E-02  -4.81482100E-04
   45.2500000       1.12214504E-04   1.04544619E-02   1.70843117E-03
   45.5000000       1.03356506E-04   9.48500354E-03   3.65940086E-03
   45.7500000       9.51617258E-05   8.18615966E-03   5.30551793E-03
   46.0000000       8.75854676E-05   6.63661631E-03   6.59854384E-03
   46.2500000       8.05809395E-05   4.92063584E-03   7.50788115E-03
   46.5000000       7.41085241E-05   3.12396954E-03   8.02180450E-03
   46.7500000       6.81304809E-05   1.33025553E-03   8.14621989E-03
   47.0000000       6.26131805E-05  -3.82547092E-04   7.90359639E-03
   47.2500000       5.75239283E-05  -1.94511819E-03   7.33078737E-03
   47.5000000       5.28342462E-05  -3.29922186E-03   6.47683442E-03
   47.7500000       4.85153796E-05  -4.39990917E-03   5.39964624E-03
   48.0000000       4.45417172E-05  -5.21612167E-03   4.16338677E-03
   48.2500000       4.08912201E-05  -5.73184993E-03   2.83498038E-03
   48.5000000       3.75380732E-05  -5.94517495E-03   1.48086634E-03
   48.7500000       3.44640721E-05  -5.86830452E-03   1.64543017E-04
   49.0000000       3.16458754E-05  -5.52532263E-03  -1.05673494E-03
   49.2500000       2.90670123E-05  -4.95122140E-03  -2.13363906E-03
   49.5000000       2.67086161E-05  -4.18899208E-03  -3.02670826E-03
   49.7500000       2.45537994E-05  -3.28746904E-03  -3.70760635E-03
   50.0000000       2.25859767E-05  -2.29871925E-03  -4.15955111E-03
   50.2500000       2.07900484E-05  -1.27538911E-03  -4.37760539E-03
   50.5000000       1.91510662E-05  -2.68222910E-04  -4.36796527E-03
   50.7500000       1.76562789E-05   6.75930351E-04  -4.14721575E-03
   51.0000000       1.62920969E-05   1.51642074E-03  -3.74066341E-03
   51.2500000       1.50459600E-05   2.22012261E-03  -3.18072550E-03
   51.5000000       1.39078247E-05   2.76246527E-03  -2.50531663E-03
   51.7500000       1.28667216E-05   3.12830973E-03  -1.75510696E-03
   52.0000000       1.19132083E-05   3.31187923E-03  -9.71938483E-04
   52.2500000       1.10383889E-05   3.31659243E-03  -1.96480338E-04
   52.5000000       1.02348858E-05   3.15441191E-03   5.33452840E-04
   52.7500000       9.49493551E-06   2.84453714E-03   1.18471228E-03
   53.0000000       8.81298001E-06   2.41236528E-03   1.73016579E-03
   53.2500000       8.18295757E-06   1.88759482E-03   2.14940519E-03
   53.5000000       7.60058947E-06   1.30286685E-03   2.42963526E-03
   53.7500000       7.06125365E-06   6.91540656E-04   2.56574061E-03
   54.0000000       6.56130624E-06   8.65185284E-05   2.56004324E-03
   54.2500000       6.09771769E-06  -4.81785974E-04   2.42190016E-03
   54.5000000       5.66807194E-06  -9.86337429E-04   2.16684351E-03
   54.7500000       5.26988970E-06  -1.40527904E-03   1.81523571E-03
   55.0000000       4.90160483E-06  -1.72213721E-03   1.39134773E-03
   55.2500000       4.56107364E-06  -1.92664238E-03   9.21478495E-04
   55.5000000       4.24703057E-06  -2.01481208E-03   4.33085253E-04
   55.7500000       3.95831739E-06  -1.98899349E-03  -4.71407257E-05
   56.0000000       3.69336885E-06  -1.85712241E-03  -4.94433974E-04
   56.2500000       3.45115473E-06  -1.63231161E-03  -8.86968686E-04
   56.5000000       3.23069253E-06  -1.33176195E-03  -1.20710512E-03
   56.7500000       3.02998728E-06  -9.75543051E-04  -1.44163205E-03
   57.0000000       2.84875432E-06  -5.86027687E-04  -1.58282206E-03
   57.2500000       2.68523104E-06  -1.85560435E-04  -1.62812718E-03
   57.5000000       2.53796611E-06   2.03555363E-04  -1.58004160E-03
   57.7500000       2.40564714E-06   5.61075518E-04  -1.44597411E-03
   58.0000000       2.28684030E-06   8.69149750E-04  -1.23750523E-03
   58.2500000       2.18004379E-06   1.11352035E-03  -9.69595916E-04
   58.5000000       2.08418101E-06   1.28415588E-03  -6.59639831E-04
   58.7500000       1.99751253E-06   1.37509475E-03  -3.26537789E-04
   59.0000000       1.91825552E-06   1.38496910E-03   1.07686501E-05
   59.2500000       1.84583598E-06   1.31707161E-03   3.33404052E-04
   59.5000000       1.77886386E-06   1.17846020E-03   6.24576176E-04
   59.7500000       1.71609474E-06   9.79659846E-04   8.69690382E-04
   60.0000000       1.65675533E-06   7.34195753E-04   1.05721899E-03
   60.2500000       1.60008653E-06   4.57279501E-04   1.17939897E-03
   60.5000000       1.54523082E-06   1.65317382E-04   1.23203127E-03
   60.7500000       1.49180823E-06  -1.25039660E-04   1.21497875E-03
   61.0000000       1.43933858E-06  -3.97983909E-04   1.13178953E-03
   61.2500000       1.38795008E-06  -6.39222562E-04   9.89618362E-04
   61.5000000       1.33744049E-06  -8.36748630E-04   7.98305788E-04
   61.7500000       1.28791748E-06  -9.81242862E-04   5.70157717E-04
   62.0000000       1.23950031E-06  -1.06667459E-03   3.18913371E-04
   62.2500000       1.19221420E-06  -1.09028269E-03   5.91419230E-05
   62.5000000       1.14660247E-06  -1.05294399E-03  -1.94708467E-04
   62.7500000       1.10284736E-06  -9.58667370E-04  -4.28724015E-04
   63.0000000       1.06149491E-06  -8.14411498E-04  -6.31053816E-04
   63.2500000       1.02220065E-06  -6.29625283E-04  -7.91057828E-04
   63.5000000       9.86017994E-07  -4.15819057E-04  -9.01727472E-04
   63.7500000       9.52671201E-07  -1.85376150E-04  -9.58283315E-04
   64.0000000       9.22579659E-07   4.84978373E-05  -9.59284953E-04
   64.2500000       8.95645258E-07   2.72924779E-04  -9.06177331E-04
   64.5000000       8.71712871E-07   4.75869922E-04  -8.03281204E-04
   64.7500000       8.50890217E-07   6.46751549E-04  -6.57725381E-04
   65.0000000       8.32600733E-07   7.77036184E-04  -4.78346657E-04
   65.2500000       8.17497892E-07   8.60961154E-04  -2.76122708E-04
   65.5000000       8.04306239E-07   8.94653786E-04  -6.24569948E-05
   65.7500000       7.93271283E-07   8.77814833E-04   1.50706299E-04
   66.0000000       7.83613416E-07   8.12345650E-04   3.51721421E-04
   66.2500000       7.75083890E-07   7.03064492E-04   5.29890764E-04
   66.5000000       7.67190897E-07   5.56849758E-04   6.76098512E-04
   66.7500000       7.59789145E-07   3.82554834E-04   7.83224707E-04
   67.0000000       7.52100675E-07   1.90005696E-04   8.46166979E-04
   67.2500000       7.44020213E-07  -9.82946949E-06   8.62510060E-04
   67.5000000       7.35217043E-07  -2.06034339E-04   8.32326186E-04
   67.7500000       7.25578332E-07  -3.87957814E-04   7.58331793E-04
   68.0000000       7.14657517E-07  -5.46270981E-04   6.45170920E-04
   68.2500000       7.02632406E-07  -6.72883703E-04   4.99859860E-04
   68.5000000       6.89391982E-07  -7.61481584E-04   3.30964860E-04
   68.7500000       6.75159754E-07  -8.08260520E-04   1.47900952E-04
   69.0000000       6.59955504E-07  -8.11426551E-04  -3.92733491E-05
   69.2500000       6.44039119E-07  -7.71705643E-04  -2.20248709E-04
   69.5000000       6.27499276E-07  -6.91883382E-04  -3.85741674E-04
   69.7500000       6.10634459E-07  -5.77153929E-04  -5.26809075E-04
   70.0000000       5.93891286E-07  -4.34259069E-04  -6.36639888E-04
   70.2500000       5.77401067E-07  -2.71414989E-04  -7.09742890E-04
   70.5000000       5.61463480E-07  -9.77483578E-05  -7.42905599E-04
   70.7500000       5.46513490E-07   7.73743959E-05  -7.35205249E-04
   71.0000000       5.32602030E-07   2.44425057E-04  -6.87647029E-04
   71.2500000       5.20218748E-07   3.94809089E-04  -6.03609544E-04
   71.5000000       5.09047084E-07   5.20398433E-04  -4.88090678E-04
   71.7500000       4.99472378E-07   6.15238445E-04  -3.47784488E-04
   72.0000000       4.91516062E-07   6.74632378E-04  -1.90754363E-04
   72.2500000       4.85039209E-07   6.95979863E-04  -2.55195191E-05
   72.5000000       4.80055860E-07   6.78784680E-04   1.38950418E-04
   72.7500000       4.76357087E-07   6.24515756E-04   2.93831865E-04
   73.0000000       4.73656655E-07   5.36622712E-04   4.30920743E-04
   73.2500000       4.71931088E-07   4.20409604E-04   5.43311005E-04
   73.5000000       4.70956223E-07   2.82444933E-04   6.25444693E-04
   73.7500000       4.70227178E-07   1.30374683E-04   6.73223287E-04
   74.0000000       4.69496285E-07  -2.76037026E-05   6.84641767E-04
   74.2500000       4.68805212E-07  -1.83052209E-04   6.59770507E-04
   74.5000000       4.67654587E-07  -3.27507558E-04   6.00327738E-04
   74.7500000       4.65807375E-07  -4.53743269E-04   5.09827863E-04
   75.0000000       4.63113281E-07  -5.55134378E-04   3.93623079E-04

Perhaps it's better to just see things graphically. gnuplot is a very good, universally available plotting program, and with just a few keystrokes we can generate good-looking graphs.

In [7]:
%%file Vo.dat
-75	0
-1	0
-1	1
1	1
1	0
75	0
Writing Vo.dat
In [8]:
%%bash
gnuplot --persist
plot 'Vo.dat' with lines
In [9]:
%%bash
gnuplot --persist
set style data line
plot [-80:80] [-0.1:1.1] 'Vo.dat' title "energy barrier V",'<./packet 40' using 1:4 title "Im packet at time = 40"

Another way of using gnuplot is to load it right into the Jupyter notebook itself.

In [10]:
%load_ext gnuplot_kernel
In [11]:
%gnuplot inline pngcairo font "Arial,16" size 800,600
In [23]:
%%gnuplot
set style data line
set yrange [-0.1:1.1]
set xrange [-60:60]
set xlabel "x/a"
set ylabel "Psi^2"
set title "Scattering of a Gaussian wave packet (41 waves), E=V/4"
plot \
  'Vo.dat' t "V(x)",\
  '<./packet  0' title "t=0",\
  '<./packet 20' title "t=20",\
  '<./packet 50' title "t=50"
set style data line
set yrange [-0.1:1.1]
set xrange [-60:60]
set xlabel "x/a"
set ylabel "Psi^2"
set title "Scattering of a Gaussian wave packet (41 waves), E=V/4"
set output '/tmp/gnuplot-inline-1517961252.4544995.585617581060.png'
plot    'Vo.dat' t "V(x)",   '<./packet  0' title "t=0",   '<./packet 20' title "t=20",   '<./packet 50' title "t=50"
gnuplot> 
gnuplot> 
unset output

Now we are ready to add the same computational functionality to the skeleton C program we introduced last time.

We already know how to compile (using a Makefile) and start a program, and how to pass command-line parameters to it at runtime.

In [13]:
%%file Makefile

OBJ = packet.o

# these are pre-defined, but we can also change the default behaviour, if we have multiple compilers
CC = cc
#CC = icc
CFLAGS = -O -DVERSION="\"`date '+%Y-%b-%d-%H:%M'`\"" 
LIBS = -lm

# special abbreviations: $@ is the file to be made; $? is the changed dependent(s).

packet: $(OBJ)
	$(CC) $(CFLAGS) -o $@ $(OBJ) $(LIBS)

# this is the default implicit rule to convert xxx.c into xxx.o 
.c.o:
	$(CC) $(CFLAGS) -c $<

clean:
	rm -f packet *.o *~ core
Writing Makefile
In [14]:
%%file packet.c
/*
 * packet.c
 *   packet - generate a Gaussian wavepacket impacting on an energy barrier.
 *
 * Completed: January.2018 (c) E.Sternin
 * Revisions: 
 *
 */

#ifndef VERSION           /* date-derived in Makefile */
#define VERSION "2018.01" /* default, that's when we first wrote the program */
#endif

#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <math.h>
#include <complex.h>

#define MAX_STR 256
#define NPTS    601       /* these cannot be changed via command-line switches */
#define X_MIN   -75.0
#define X_MAX    75.0
#define NS        4       /* +/- NS * sigma is the extent of the Gaussian packet we calculate */

/* Global variables */

static char whoami[MAX_STR] ;        /* argv[0] will end up here */
static int  verbosity;               /* Show detailed information */
static char options[] = "Vvhn:t:U:p:s:x:";   /* command-line options, : means takes a value */
static char help_msg[] = "\
%s [<options>]\n\
\t-V        report version information\n\
\t-v        increase verbosity, may combine several\n\
\t-n #      number of waves making up the packet\n\
\t-t #      time since the beginning\n\
\t-U #      U[=16] height of the barrier U(x), in epsilon\n\
\t-p #      p0[=0.5] is the mean momentum of the packet\n\
\t-s #      s[=0.1] is the width of the packet in p-space\n\
\t-x #      x0[=-25] is the initial position of the packet\n\
\t-h        help message\n\
\n\
e.g.\tpacket -v -n 20 -t 20 -p 0.5 -s 0.1\n" 
;

/*************************************************service routines************/

void __attribute__((noreturn)) die(char *msg, ...) {
  va_list args;
  va_start(args, msg);
  vfprintf(stderr, msg, args);
  fputc('\n', stderr);
  exit(1);
  }

/************************************************************** main *********/
int main(int argc, char **argv) {
  int    i,j,n;
  double t,U,p0,s,x0,dx,x[NPTS],A0,p_step,p_now,E;
  double complex Psi[NPTS];

  double Pi = 2.*acos(0.);

/*
 *  default values, may get changed by the command-line options
 */
  verbosity = 0;    /* 0 = quiet by default, 1 = info, 2 = debug */	
  n = 20;           /* default to 41=(2*n+1) waves in the packet */
  p0 = 0.5;         /* default to 0.5 in p-space */
  s = 0.1;          /* keep to less than 1/4 of p0, so that (p0+/-4s) > 0, all */
                    /* constituent waves travel to the right initially */
  U = 1;            /* barrier height, in epsilon */
  x0 = -25;         /* initial position of the packet */
  t = 0;

  strncpy(whoami, argv[0], MAX_STR);

  while ((i = getopt(argc, argv, options)) != -1)
    switch (i) {
      case 'V':
         printf(" %s: version %s\n",whoami,VERSION);
         break;
      case 'v':
         verbosity++;
         if (verbosity > 1) printf(" %s: verbosity level set to %d\n",whoami,verbosity);
         break;
      case 'h':
         die(help_msg,whoami);
         break;
      case 'n':
         if ( ((n = atoi(optarg)) > 10000) || n < 10 )
           die(" %s: -n %d is not a valid number of points (10..10000)\n",whoami,n);
         if (verbosity > 1) printf(" %s: Number of points = %d\n",whoami,n);
         break;
      case 't':
         if ( ((t = atof(optarg)) < 0) )
           die(" %s: -t %d is not a valid time\n",whoami,t);
         if (verbosity > 1) printf(" %s: time = %f\n",whoami,t);
         break;
      case 'U':
         U = atof(optarg);
         if (verbosity > 1) printf(" %s: barrier height is = %f epsilon\n",whoami,U);
         break;
      case 'p':
         p0 = atof(optarg);
         if (verbosity > 1) printf(" %s: mean momentum of packet is p0 = %f\n",whoami,p0);
         break;
      case 's':
         if ( (s = atof(optarg)) < 1e-2 )
           die(" %s: -s %f is not a valid packet width, >= 1e-2\n",whoami,s);
         if (verbosity > 1) printf(" %s: width of packet, s = %f\n",whoami,s);
         break;
      case 'x':
         if ( ((x0 = atof(optarg)) < X_MIN) || x0 > 0)
           die(" %s: -x %f is not a valid packet position, %d..0\n",whoami,x0,X_MIN);
         if (verbosity > 1) printf(" %s: initial position of packet, x0 = %f\n",whoami,x0);
         break;
      default:
         if (verbosity > 0) die(" try %s -h\n",whoami);	/* otherwise, die quietly */
         return 0;
      }

/*
 * when we get here, we parsed all user input, and are ready to calculate things
 */

  if ((p0-NS*s) <= 0) /* is the smallest value of p still positive ? */
     die(" %s: some waves in the packet p=%f+/-%d*%f are not travelling to the right at t=0\n",whoami,p0,NS,s);

  if ((p0+NS*s) >= sqrt(U)) /* is the largest value of p still below U of the barrier ? */
     die(" %s: some waves in the packet p=%f+/-%d*%f have energy >= U=%f\n",whoami,p0,NS,s,U);

  dx = (X_MAX-X_MIN)/(double)(NPTS-1);
  for (j=0; j < NPTS; j++) { 
     x[j]=X_MIN + j*dx;
     Psi[j] = (double complex)0.0;  /* start with a blank Psi(x) */
     }

  A0 = 1./sqrt( sqrt(2.*Pi) * s ) /(double)(2*n+1);
  p_step = NS*s/(double)n;    /* calculate over p0 +/- NS*sigma_p */

  double complex A,k1,k2,ce4,C1,b,c,d,f;

  /* add (complex) waves to Psi[j]; each wave is calculated at x[j]*/
  for (i=-n; i < n; i++) {
    p_now = p0 + (double)i * p_step;
    E = p_now*p_now;
    A = A0*exp(-pow((p_now-p0)/(2*s),2))*cexp(-I*(p_now*x0 + E*t));

    k1 = csqrt((double complex)(E));
    k2 = csqrt((double complex)(E-U));
    ce4 = cexp(4*I*k2);
    C1 = -2*k1*k2-k1*k1-2*k1*k2*ce4+ce4*(k1*k1+k2*k2)-k2*k2;

    b = (k1*k1-k2*k2)/ C1 *(cexp(2*I*(2*k2-k1))-cexp(-2*I*k1));
    c = -2*k1*(k2+k1)/ C1 *cexp(I*(k2-k1));
    d = (-2*k2+2*k1)*k1/ C1 *cexp(I*(-k1+3*k2));
    f = -4*k1*k2/ C1 *cexp(2*I*(k2-k1));

    for (j=0; j < NPTS; j++) {
       if (x[j] < -1.) {         /* to the left of the barrier, x < -1 */
          Psi[j] += A*( cexp(I*k1*x[j]) + b*cexp(-I*k1*x[j]) );
          }
       else if (x[j] > 1.) {     /* to the right of the barrier, x > 1 */
          Psi[j] += A*f*cexp(I*k1*x[j]);
          }
       else {                    /* inside the barrier, -1 < x < 1 */
          Psi[j] += A*( c*cexp(I*k2*x[j]) + d*cexp(-I*k2*x[j]) );
          }
       }
  }

  /* output the total Psi as is, without normalization */
  for (j=0; j < NPTS; j++) {
     if (verbosity > 0) printf("\t%f\t%f\t%f\t%f\n",x[j],cabs(Psi[j]),creal(Psi[j]),cimag(Psi[j]));
     }

  return 0;
  }
Writing packet.c
In [24]:
%%bash
make clean; make
./packet -x -40 -t 0
rm -f packet *.o *~ core
cc -O -DVERSION="\"`date '+%Y-%b-%d-%H:%M'`\""  -c packet.c
cc -O -DVERSION="\"`date '+%Y-%b-%d-%H:%M'`\""  -o packet packet.o -lm

The homework for today was to convert the evolution expressions into dimensionless form, ready to be implemented in code. Today's session is to try to do just that.

In [27]:
%%gnuplot
set style data line
set xrange [-75:75]
set yrange [-0.05:1.8]
set xlabel "x/a"
set ylabel "Psi(x)"
set title "Scattering of a Gaussian wave packet (51 waves), E=p_0^2=U/4"
plot \
  'Vo.dat' t "U(x)",\
  '<./packet -v -n 25 -x -40 -p .5 -s 0.07 -t 0' title "t=0",\
  '<./packet -v -n 25 -x -40 -p .5 -s 0.07 -t 40' title "t=40",\
  '<./packet -v -n 25 -x -40 -p .5 -s 0.07 -t 70' title "t=70"
set style data line
set xrange [-75:75]
set yrange [-0.05:1.8]
set xlabel "x/a"
set ylabel "Psi(x)"
set title "Scattering of a Gaussian wave packet (51 waves), E=p_0^2=U/4"
set output '/tmp/gnuplot-inline-1517961329.5172246.5513159542.png'
plot    'Vo.dat' t "U(x)",   '<./packet -v -n 25 -x -40 -p .5 -s 0.07 -t 0' title "t=0",   '<./packet -v -n 25 -x -40 -p .5 -s 0.07 -t 40' title "t=40",   '<./packet -v -n 25 -x -40 -p .5 -s 0.07 -t 70' title "t=70"
gnuplot> 
gnuplot> 
unset output

Notes on Fortran-to-C conversion

Fortran:

    real*4      x(601), ...
    complex*8   P(601), ...
    ...
    do j = 1,601
      x(j) = -75 + (j-1)*0.25
      P(j) = cmplx(0.,0.)
    end do

C:

...
#include <math.h>
#include <complex.h>

#define NPTS    601       /* these cannot be changed via command-line switches */
#define X_MIN   -75.0
#define X_MAX    75.0
...
  double dx,x[NPTS], ...;
  double complex Psi[NPTS];
...
  dx = (X_MAX-X_MIN)/(double)(NPTS-1);
  for (j=0; j < NPTS; j++) { 
     x[j]=X_MIN + j*dx;
     Psi[j] = (double complex)0.0;  /* start with a blank Psi(x) */
     }

Fortran:

    I = cmplx(0.,1.)    !complex i

    k1 = sqrt(cmplx(E,0.))
    K2 = sqrt(cmplx(E-1.,0.))

    C1 = -2*k1*K2-k1**2-2*k1*K2*exp(4*I*K2)+exp(4*I*K2)*k1**2
     *  -K2**2+K2**2*exp(4*I*K2)

    b = (k1**2-K2**2)/ C1 *(exp(2*I*(2*K2-k1))-exp(-2*I*k1))
    c = -2*k1*(K2+k1)/ C1 *exp(I*(K2-k1))
    d = (-2*K2+2*k1)*k1/ C1 *exp(I*(-k1+3*K2))
    f = -4*k1*K2/ C1 *exp(2*I*(K2-k1))

C:

#include <complex.h>
...
  for (i=-n; i < n; i++) {
    p_now = p0 + (double)i * p_step;
    E = p_now*p_now;
    A = A0*exp(-pow((p_now-p0)/(2*s),2))*cexp(-I*(p_now*x0 + E*t));

    k1 = csqrt((double complex)(E));
    k2 = csqrt((double complex)(E-U));
    ce4 = cexp(4*I*k2);
    C1 = -2*k1*k2-k1*k1-2*k1*k2*ce4+ce4*(k1*k1+k2*k2)-k2*k2;

    b = (k1*k1-k2*k2)/ C1 *(cexp(2*I*(2*k2-k1))-cexp(-2*I*k1));
    c = -2*k1*(k2+k1)/ C1 *cexp(I*(k2-k1));
    d = (-2*k2+2*k1)*k1/ C1 *cexp(I*(-k1+3*k2));
    f = -4*k1*k2/ C1 *cexp(2*I*(k2-k1));

    ...
  }

Fortran:

    do j = 1,296        ! left of the barrier, x < -1
      P(j) = P(j) + A * ( exp(I*k1*x(j)) + b*exp(-I*k1*x(j)) )
    end do
    do j=297,304        ! inside the barrier, -1 < x < 1
      P(j) = P(j) + A * ( c*exp(I*K2*x(j)) + d*exp(-I*K2*x(j)) )
    end do
    do j=305,601        ! to the right of the barrier, x > 1
      P(j) = P(j) + A * f*exp(I*k1*x(j))
    end do

C:

    for (j=0; j < NPTS; j++) {
       if (x[j] < -1.) {         /* to the left of the barrier, x < -1 */
          Psi[j] += A*( cexp(I*k1*x[j]) + b*cexp(-I*k1*x[j]) );
          }
       else if (x[j] > 1.) {     /* to the right of the barrier, x > 1 */
          Psi[j] += A*f*cexp(I*k1*x[j]);
          }
       else {                    /* inside the barrier, -1 < x < 1 */
          Psi[j] += A*( c*cexp(I*k2*x[j]) + d*cexp(-I*k2*x[j]) );
          }
       }