can be issued as a commandline parameter |'
+ Write(6,*) ' ---------------------------------------------------------------------'
+
+! files :
+
+ inquire( file='\simulations\counter.his', exist = unique )
+ IF ( unique ) THEN
+ open(9,file='\simulations\counter.his',status='old',err=994)
+ read(9,*) ifile ! initialize outputfile counter
+ ELSE
+ ifile = 0
+ END IF
+
+ IF ( iargc() .GT. 0 ) THEN
+ call getarg(1, calculation)
+ Write(6,*) ' Calculation taken from commandline > ',calculation
+ ELSE
+200 write(6,201)
+201 format(' '/' Give name of the calculation > ', \)
+ read(5,'(a60)') calculation
+ END IF
+
+ l_calc = index( calculation, ' ') - 1
+
+ IF ( l_calc .GT. 0 ) THEN
+ open(1,file=calculation(1:l_calc)//'.inp',status='old',action='read',err=995 )
+ open(2,file=calculation(1:l_calc)//'.out',status='unknown',action='write',err=996)
+
+
+ END IF
+
+! initialization of randomumber generator and Pi
+
+ iseed = 1234567
+ Pi = acos( -1.0D+00 )
+
+! Read everything from the input file, one line per calculation
+
+ DO WHILE ( .NOT. Eof(1) ) ! WHILE LOOP(1) over the input file
+
+10 read(1,'(a127)',end=999) line
+
+ IF ( line(1:1) .EQ. '!' ) THEN
+ Write(2,'(a)') line
+ GOTO 10
+ END IF
+
+ read(line,*,err=998,end=999) a, moment, b_ext, d, w, concentration, &
+ & n_site, t_ini, p_ini, depth1, depth2, &
+ & fluctuationrate, dphi, dtheta
+
+ IF ( n_site .GT. 0.8 * max_muons ) n_site = 0.8 * max_muons
+
+! Estimate optimum time_step
+
+ b_abs = sqrt( sum( b_ext * b_ext ) )
+ f_c = 135.5 * b_abs + 0.3 * moment * comcentration / (a*a*a)
+
+ time_step = min( 0.01, 0.1 / f_c )
+
+! Initialize randomnumber generator "randomly"
+
+ call date_and_time( dddd, tttt, zone, dt )
+ DO i = 1, dt(8) ! number milliseconds on the clock
+ dummy = ran(iseed)
+ END DO
+!
+
+ write(2,100) calculation(1:73),(dt(j),j=1,3),(dt(j),j=5,8)
+100 format(' '/' ',73('-')/' ',a73/' ',73('-')/ &
+& ' Calculation started ',i5,'-',i2,'-',i2, &
+& ' at ',2(i2,':'),i2,'.',i3/' ',73('-')/' ')
+ write(2,'(a,f8.3)') ' lattice parameter = ', a
+ write(2,'(a,f8.3)') ' magnetic moment = ', moment
+ write(2,'(a,3f8.3)') ' external field = ', b_ext
+ write(2,'(a,f8.3)') ' concentration = ', concentration
+ write(2,'(a,2f8.3)') ' init.muon theta,phi = ', t_ini, p_ini
+ write(2,'(a,f8.3)') ' fluctuationrate = ', fluctuationrate
+ write(2,'(a,2f8.3)') ' fluctuation amp. = ', dphi, dtheta
+
+ emu(1) = sin(t_ini*Pi/180.0) * cos(p_ini*Pi/180.0)
+ emu(2) = sin(t_ini*Pi/180.0) * sin(p_ini*Pi/180.0)
+ emu(3) = cos(t_ini*Pi/180.0)
+
+ DO j = 1, max_muons
+ m(j).dir = emu
+ END DO
+
+ exp_time = 0.0D+00
+ write_time = 0.0D+00
+
+! update file index for
+
+ ifile = ifile + 1 ! increase outputfile number
+
+ IF ( unique ) THEN
+ rewind(9)
+ write(9,*) ifile ! store for next program
+ END IF
+
+ write(file_index,'(''_'',i3)') ifile ! generate file_name
+
+ DO j = 2, 4
+ IF (file_index(j:j) .EQ. ' ' ) file_index(j:j) = '0'
+ END DO
+
+ open(3,file=calculation(1:l_calc)//file_index//'.g_t', &
+ & status='unknown',action='write', err=997 )
+
+! output time=0 asymmetries
+
+ write(3,'(4F19.6)' ) exp_time, emu
+ write_time = write_time + time_resolution
+
+! Initialize randomnumber generator "randomly"
+
+ call date_and_time( dddd, tttt, zone, dt )
+
+ DO i = 1, dt(8) ! number milliseconds on the clock
+ dummy = ran(iseed)
+ END DO
+
+
+! make lattice, spinglass, choose muon sites and calculates "interaction matrix"
+
+ CALL lattice( iseed, d, w, a, concentration, n_spin, nd, nw, s, &
+ & n_site, depth1, depth2, m )
+
+
+ write(2,'(a,f8.3)') ' thickness = ', d
+ write(2,'(a,f8.3)') ' width = ', w
+ write(2,'(a,i10)') ' number of spins = ', n_spin
+ write(2,'(a,i10)') ' number of muons = ', n_site
+ write(2,'(a,2f8.3)') ' muons penetrate betw. ', depth1, depth2
+ write(2,'(a,a)') ' Output will be written on ', &
+ & calculation(1:l_calc)//file_index//'.g_t'
+
+ write(6,'(a,f8.3)') ' lattice parameter = ', a
+ write(6,'(a,f8.3)') ' magnetic moment = ', moment
+ write(6,'(a,3f8.3)') ' external field = ', b_ext
+ write(6,'(a,f8.3)') ' concentration = ', concentration
+ write(6,'(a,2f8.3)') ' init.muon theta,phi = ', t_ini, p_ini
+ write(6,'(a,f8.3)') ' fluctuationrate = ', fluctuationrate
+ write(6,'(a,2f8.3)') ' fluctuation amp. = ', dphi, dtheta
+ write(6,'(a,f8.3)') ' thickness = ', d
+ write(6,'(a,f8.3)') ' width = ', w
+ write(6,'(a,i10)') ' number of spins = ', n_spin
+ write(6,'(a,i10)') ' number of muons = ', n_site
+ write(6,'(a,2f8.3)') ' muons penetrate betw. ', depth1, depth2
+ write(6,'(a,a)') ' Output will be written on ', &
+ & calculation(1:l_calc)//file_index//'.g_t'
+
+
+
+! The fluctuations are incorporated as follows:
+
+! for rates larger then 1/time_resolution MHz,
+! a timestep tau is choosen from a
+! log distrubution ( tau = - ln(random) / fluctuationrate )
+! The muon then rotates for tau microseconds, after all spins
+! are rotated over an angle between - dtetha en dtheta and - dphi and dphi.
+! This process is repeated until the total time is n_time_steps*time_resolution
+! microsecods or more.
+! Output of the muon position is about every time_resolution microsecond.
+
+! for rates smaller then 1/time_resolution MHz
+! the muons rotate 1000 times for time_resolution microsecond, after each rotation
+! a fraction (= fluctuationrate *time_resolution) of the magnetic ions are rotated
+! over an angle between - dtetha en dtheta and - dphi and dphi.
+
+! After each fluctuation the fields at the muonsites are recalculated.
+
+! "deporization" functions are calculated for
+! left-right, up-down and forward-backward detectors,
+! being the x-, y- and z-components of the muon spin vector.
+
+! For arbitrary direction one has to take the scalar product of
+! that specific direction with the results produced by this program
+
+
+ IF ( fluctuationrate .GT. 1.0/time_resolution ) THEN ! Rapid fluctuations
+
+ fraction = 1.0
+ ! Start of WHILE loop(2) over exp_time
+
+ DO WHILE ( exp_time .LT. time_resolution * float(n_time_steps) )
+
+ tau = - log( ran(iseed) ) / fluctuationrate
+
+ IF ( exp_time + tau .GT. time_resolution * float(n_time_steps) ) &
+ & tau = time_resolution * float(n_time_steps) - exp_time
+
+! take at least time_resolution microsec. steps, even if tau is larger
+
+ n_steps = floor( tau / time_resolution ) + 1
+ step = tau / float( n_steps )
+
+ CALL fields( a, moment, b_ext, s, n_site, m)
+
+ DO i_step = 1, n_steps
+ exp_time = exp_time + step
+ call muonrotation( n_site, m, step )
+ g_t = 0.0
+ DO j = 1, n_site
+ DO k = 1, 3
+ g_t(k) = g_t(k) + m(j).dir(k)
+ END DO
+ END DO
+
+ g_t = g_t / float(n_site)
+
+ IF ( exp_time .GT. write_time ) THEN
+ write(3,'(4F19.6)' ) exp_time, g_t
+ write_time = exp_time + time_resolution
+ END IF
+
+ END DO
+
+! after tau, change spin directions and repeat the above.
+! however, stop when 10 microsec has been reached.
+
+ CALL fluctuation( iseed, n_spin, s, dtheta, dphi, fraction )
+
+ END DO ! END of WHILE loop(2)
+
+ ELSE ! fluctuationrate < 1/time_resolution
+
+ fraction = fluctuationrate * time_resolution
+ step = time_resolution
+ n_steps = n_time_steps
+
+ DO i_step = 1, n_steps
+ exp_time = exp_time + step
+ CALL fields( a, moment, b_ext, s, n_site, m)
+ CALL muonrotation( n_site, m, step )
+
+ g_t = 0.0
+
+ DO j = 1, n_site
+ DO k = 1, 3
+ g_t(k) = g_t(k) + m(j).dir(k)
+ END DO
+ END DO
+
+ g_t = g_t / float(n_site)
+
+ write(3,'(4F19.6)' ) exp_time, g_t
+ CALL fluctuation( iseed, n_spin, s, dtheta, dphi, fraction )
+
+ END DO
+
+ END IF
+
+ END DO ! END of WHILE loop(1)
+
+ STOP ' Program DYNAMICS stopped where it should stop '
+
+994 STOP ' FATAL: Cannot open counter.his '
+
+995 STOP ' FATAL: Cannot open input file '
+
+996 STOP ' FATAL: Cannot open output file '
+
+997 Write(2,*) ' Cannot open g_t file '
+ STOP ' FATAL: Cannot open g_t file '
+
+998 Write(2,*) ' Error in input file '
+ STOP ' FATAL: Due to error in input file '
+
+999 Write(2,*) ' End of input-file '
+ STOP ' STOP End of input file '
+
+ end program dynamics
+
+!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! LATTICE calculates the actual dimensions of the sample, places magnetic spins
+! randomly according to concentration, gives the spins a random direction
+! in space. It also generates a table of muonsites.
+
+!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+ Subroutine lattice( iseed, d, w, a, concentration, n_spin, nd, nw, s, &
+ & n_site, depth1, depth2, m )
+
+! Structure to store the position (as lattice site-indexes)
+! and the direction-cosines of each spin and muon.
+
+ implicit none
+
+ integer*4,parameter::max_spins = 50000, & ! maximum number of magnetic moments
+ & max_muons = 10000, & ! maximum number of muons
+ & max_nn = 500 ! maximum number of nearest neighbours
+
+
+ structure /spin/
+ integer*4 x,y,z
+ real*8 theta,phi,dir(3)
+ end structure
+
+ structure /muon/
+ integer*4 x,y,z, ns, s(max_nn)
+ real*8 dir(3), r(3,max_nn), r_2(max_nn), r_5(max_nn), omega(3)
+ end structure
+
+ real*8 d, w, a, concentration, c, depth1, depth2, fraction, radiussq
+ real*8 Pi, r(3), r_2, r_3, r_5, help
+ integer*4 iseed, nd, nw, nat, n_spin, n_site, nd1, nd2
+ integer*4 i, j, k, l, hw, kw, ns
+ record /spin/ s(*)
+ record /muon/ m(*)
+
+ Pi = acos( -1.0D+00 )
+ c = concentration / 100.0
+
+! Calculate the 'rounded' number of spins for a lattice m*m*n for
+! the given concentration.
+! n is the number of atoms (half unitcells) perpendicular
+! to the layer (== z-direction).
+! m is the size of the layer in the x- and y-direction
+
+ nd = floor(2.0 * d / a ) + 2
+ nw = floor(2.0 * w / a ) + 2
+ nat = nd * nw * nw / 2
+ n_spin = floor( nat * c )
+ d = float(nd) * a / 2.0
+ w = float(nw) * a / 2.0
+ hw = nw / 2
+
+ IF ( c .GT. 0.0 ) THEN
+ radiussq = (( 1.6 * float(max_nn) / c ) / ( 4.0 * Pi / 3.0 ))**(2.0/3.0)
+ ELSE
+ radiussq = 1D+10
+ END IF
+
+ write(2,*) ' radius = ', sqrt( radiussq )
+
+ nd1 = floor( 2.0 * depth1 / a )
+ nd2 = floor( 2.0 * depth2 / a )
+ IF ( mod( nd1 , 2 ) .EQ. 0 ) nd1 = nd1 + 1 ! nd1 should be odd
+ IF ( nd2 .LT. nd1 + 1 ) nd2 = nd1 + 1
+ depth1 = float(nd1) * a / 2.0
+ depth2 = float(nd2) * a / 2.0
+ fraction = float(n_site) / (float((nd2-nd1)*nw*nw) / 8.0)
+
+
+! Place the spins randomly on the fcc-lattice
+! Run over a whole simple cubic lattice in steps
+! of half of the fcc-unitcell.
+! Then take care of the fcc-structure and
+! decide whether or not to place a spin.
+
+ n_spin = 0
+
+ DO j = 0, nw-1
+ DO k = 0, nw-1
+ DO l = 0, nd-1
+ IF ( mod(j+k+l,2) .EQ. 0 ) THEN ! This takes care of the fcc structure.
+ IF ( ran(iseed) .LT. c ) THEN ! Takes care of concentration
+ n_spin = n_spin + 1
+ IF ( n_spin .GT. max_spins ) STOP ' Stopped because number of spin too large '
+ s(n_spin).x = j
+ s(n_spin).y = k
+ s(n_spin).z = l
+
+! Give the spin an arbitrary direction
+
+ s(n_spin).theta = Pi * ran(iseed)
+ s(n_spin).phi = (Pi+Pi) * ran(iseed)
+ s(n_spin).dir(1) = sin(s(n_spin).theta) * cos(s(n_spin).phi)
+ s(n_spin).dir(2) = sin(s(n_spin).theta) * sin(s(n_spin).phi)
+ s(n_spin).dir(3) = cos(s(n_spin).theta)
+
+ END IF
+ END IF
+
+
+ END DO
+ END DO
+ END DO
+
+! determine positions of the muons
+
+ n_site = 0
+
+ DO j = 1, nw-1, 2
+ DO k = 1, nw-1, 2
+ DO l = nd1, nd2, 2
+
+ IF ( ran(iseed) .LT. fraction ) THEN
+ n_site = n_site + 1
+ m(n_site).x = j
+ m(n_site).y = k
+ m(n_site).x = l
+ ns = 0
+
+ DO i = 1, n_spin
+ kw = j - s(i).x
+ IF ( kw .LT. -hw ) kw = kw + nw ! periodic boundary condition
+ IF ( kw .GT. hw ) kw = kw - nw ! periodic boundary condition
+ r(1) = dble(float(kw))
+ kw = k - s(i).y
+ IF ( kw .LT. -hw ) kw = kw + nw ! periodic boundary condition
+ IF ( kw .GT. hw ) kw = kw - nw ! periodic boundary condition
+ r(2) = dble(float(kw))
+ r(3) = dble(float(l-s(i).z)) ! NO periodic boundary condition
+ r_2 = sum( r * r )
+
+ IF ( r_2 .LE. radiussq ) THEN ! skip calculation if distance
+ ! is too large
+ help = sqrt( r_2 )
+ r_3 = r_2 * help
+ r_5 = r_2 * r_3
+ ns = ns + 1
+ IF (ns .GT. max_nn) STOP ' Stopped because NS becomes too large '
+ m(n_site).s(ns) = i
+ m(n_site).r(1,ns) = r(1)
+ m(n_site).r(2,ns) = r(2)
+ m(n_site).r(3,ns) = r(3)
+ m(n_site).r_2(ns) = r_2
+ m(n_site).r_5(ns) = r_5
+ END IF
+
+ END DO
+
+ m(n_site).ns = ns
+
+ END IF
+
+ END DO
+ END DO
+ END DO
+
+
+
+ RETURN
+ END
+
+!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! FIELDS calculates all internal fields at the muonsites and
+! adds the external field
+
+!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+ Subroutine fields( a, moment, b_ext, s, n_site, m)
+
+ implicit none
+
+ integer*4,parameter::max_spins = 50000, & ! maximum number of magnetic moments
+ & max_muons = 10000, & ! maximum number of muons
+ & max_nn = 500 ! maximum number of nearest neighbours
+
+
+ structure /spin/
+ integer*4 x,y,z
+ real*8 theta,phi,dir(3)
+ end structure
+
+ structure /muon/
+ integer*4 x,y,z, ns, s(max_nn)
+ real*8 dir(3), r(3,max_nn), r_2(max_nn), r_5(max_nn), omega(3)
+ end structure
+
+ real*8 Pi, Gyro, p_r, a, factor
+ real*8 b(3), b_ext(3), moment
+ integer*4 j, k, l, n_site
+ record /spin/ s(*)
+ record /muon/ m(*)
+
+ Pi = acos(-1D+00)
+ Gyro = (Pi+Pi) * 135.54 ! gyro-magnetic ratio of muon [tesla^-1 s^-1]
+ factor = 1D-07 * moment * 9.2740019D-24 / ( a*a*a * 0.125D-27 )
+
+ DO j = 1, n_site
+ b = 0
+ DO k = 1, m(j).ns
+ l = m(j).s(k)
+ p_r = m(j).r(1,k) * s(l).dir(1) + &
+ & m(j).r(2,k) * s(l).dir(2) + &
+ & m(j).r(3,k) * s(l).dir(3)
+ b(1) = b(1) + (3.0D+00*p_r*m(j).r(1,k)-m(j).r_2(k)*s(l).dir(1))/m(j).r_5(k)
+ b(2) = b(2) + (3.0D+00*p_r*m(j).r(2,k)-m(j).r_2(k)*s(l).dir(2))/m(j).r_5(k)
+ b(3) = b(3) + (3.0D+00*p_r*m(j).r(3,k)-m(j).r_2(k)*s(l).dir(3))/m(j).r_5(k)
+ END DO
+ b = factor * b + b_ext
+ m(j).omega(1) = Gyro * b(1)
+ m(j).omega(2) = Gyro * b(2)
+ m(j).omega(3) = Gyro * b(3)
+ END DO
+
+ RETURN
+ END
+
+!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! FLUCTUATION changes all directions of the spins with a random amount
+! DTHETA and DPHI
+
+!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+ Subroutine fluctuation( iseed, n_spin, s, dtheta, dphi, fraction )
+
+ implicit none
+
+ integer*4,parameter::max_spins = 50000, & ! maximum number of magnetic moments
+ & max_muons = 10000, & ! maximum number of muons
+ & max_nn = 500 ! maximum number of nearest neighbours
+
+ structure /spin/
+ integer*4 x,y,z
+ real*8 theta,phi,dir(3)
+ end structure
+
+ record /spin/ s(*)
+ real*8 dtheta, dphi, dt, dp, Pi, fraction
+ integer*4 iseed, n_spin, i_spin
+
+ IF ( fraction .LT. 1.0D-06 .OR. &
+ & ( dtheta .LT. 1.0D-06 .AND. dphi .LT. 1.0D-06 ) ) RETURN
+
+ Pi = acos( -1.0D+00 )
+ dt = dtheta * Pi / 180.0 ! amplitude of the fluctuation in theta
+ dp = dphi * Pi / 180.0 ! amplitude of the fluctuation in phi
+
+ DO i_spin = 1, n_spin
+ IF ( ran(iseed) .LT. fraction ) THEN
+ s(i_spin).theta = s(i_spin).theta + 2.0 * dt * (ran(iseed)-0.5)
+ s(i_spin).phi = s(i_spin).phi + 2.0 * dp * (ran(iseed)-0.5)
+ s(i_spin).dir(1) = sin(s(i_spin).theta) * cos(s(i_spin).phi)
+ s(i_spin).dir(2) = sin(s(i_spin).theta) * sin(s(i_spin).phi)
+ s(i_spin).dir(3) = cos(s(i_spin).theta)
+ END IF
+ END DO
+
+ RETURN
+ END
+
+!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! MUONROTATION rotates all muons over the vector m.omega * step
+
+!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+ Subroutine muonrotation( n_site, m, step )
+
+ implicit none
+
+ integer*4,parameter::max_spins = 50000, & ! maximum number of magnetic moments
+ & max_muons = 10000, & ! maximum number of muons
+ & max_nn = 500 ! maximum number of nearest neighbours
+
+ structure /muon/
+ integer*4 x,y,z, ns, s(max_nn)
+ real*8 dir(3), r(3,max_nn), r_2(max_nn), r_5(max_nn), omega(3)
+ end structure
+
+ record /muon/ m(*)
+ real*8 v(3), OM(3), step
+ integer*4 j, n_site
+
+ DO j = 1, n_site
+ OM(1) = m(j).omega(1)
+ OM(2) = m(j).omega(2)
+ OM(3) = m(j).omega(3)
+ OM = OM * step
+ v(1) = m(j).dir(1)
+ v(2) = m(j).dir(2)
+ v(3) = m(j).dir(3)
+ call rotation( v, OM )
+ m(j).dir(1) = v(1)
+ m(j).dir(2) = v(2)
+ m(j).dir(3) = v(3)
+ END DO
+
+ RETURN
+ END
+
+!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! ROTATION rotates a vector V around the vector O
+
+!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+ Subroutine rotation( v, o )
+
+ implicit none
+
+ real*8 v(3), o(3), uo(3), r(3), o_abs, cc, ss
+
+ o_abs = sqrt( sum( o * o ) )
+
+ IF ( o_abs .GT. 1D-08 ) THEN
+ uo = o / o_abs
+ cc = cos( o_abs )
+ ss = sin( o_abs )
+
+ r(1) = ( cc+uo(1)*uo(1)*(1-cc) ) * v(1) + &
+ & ( -uo(3)*ss+uo(1)*uo(2)*(1-cc) ) * v(2) + &
+ & ( uo(2)*ss+uo(1)*uo(3)*(1-cc) ) * v(3)
+
+ r(2) = ( uo(3)*ss+uo(1)*uo(2)*(1-cc) ) * v(1) + &
+ & ( cc+uo(2)*uo(2)*(1-cc) ) * v(2) + &
+ & ( -uo(1)*ss+uo(2)*uo(3)*(1-cc) ) * v(3)
+
+ r(3) = ( -uo(2)*ss+uo(1)*uo(3)*(1-cc) ) * v(1) + &
+ & ( uo(1)*ss+uo(2)*uo(3)*(1-cc) ) * v(2) + &
+ & ( cc+uo(3)*uo(3)*(1-cc) ) * v(3)
+
+ v = r
+ END IF
+
+ RETURN
+ END
diff --git a/dynamics/dynamics.plg b/dynamics/dynamics.plg
new file mode 100755
index 0000000..ca5ade1
--- /dev/null
+++ b/dynamics/dynamics.plg
@@ -0,0 +1,19 @@
+
+
+
+Build Log
+
+--------------------Configuration: dynamics - Win32 Debug--------------------
+
+Command Lines
+Creating command line "link.exe kernel32.lib /nologo /subsystem:console /incremental:no /pdb:"Debug/dynamics.pdb" /debug /machine:I386 /out:"Debug/dynamics.exe" /pdbtype:sept .\Debug\dynamics.obj "
+Output Window
+Linking...
+
+
+
+Results
+dynamics.exe - 0 error(s), 0 warning(s)
+
+
+
diff --git a/index/Debug/DF60.PDB b/index/Debug/DF60.PDB
new file mode 100755
index 0000000..3247b40
Binary files /dev/null and b/index/Debug/DF60.PDB differ
diff --git a/index/Debug/index old files.obj b/index/Debug/index old files.obj
new file mode 100755
index 0000000..c879639
Binary files /dev/null and b/index/Debug/index old files.obj differ
diff --git a/index/Debug/index.pdb b/index/Debug/index.pdb
new file mode 100755
index 0000000..1c0698e
Binary files /dev/null and b/index/Debug/index.pdb differ
diff --git a/index/Release/index old files.obj b/index/Release/index old files.obj
new file mode 100755
index 0000000..95c23ed
Binary files /dev/null and b/index/Release/index old files.obj differ
diff --git a/index/index old files.f90 b/index/index old files.f90
new file mode 100755
index 0000000..b9ad32b
--- /dev/null
+++ b/index/index old files.f90
@@ -0,0 +1,60 @@
+!
+! Program to index old files
+!
+ character*4 file_index
+ character*80 calculation, dummy
+ character*48 g_t
+ character*32 his
+ integer ifile, i,j
+!
+ ifile = 0
+!
+1 CONTINUE
+200 write(6,201)
+201 format(' '/' Give name of the calculation > ', \)
+ read(5,'(a60)') calculation
+ IF ( calculation(1:1) .EQ. ' ' ) STOP ' Stopped by operator '
+ l_calc = index( calculation, ' ') - 1
+!
+ open(1,file=calculation(1:l_calc)//'.inp',status='old',action='read',err=200 )
+ open(3,file=calculation(1:l_calc)//'.g_t',status='old',action='read')
+ open(4,file=calculation(1:l_calc)//'.his',status='old',action='read')
+!
+ read(1,'(a80)') comment
+ read(4,'(a80)') dummy
+!
+300 read(1,'(a80)',end=900) dummy
+ ifile = ifile + 1 ! increase outputfile number
+ write(file_index,'(''_'',i3)') ifile ! generate file_name
+ DO j = 2, 4
+ IF (file_index(j:j) .EQ. ' ' ) file_index(j:j) = '0'
+ END DO
+!
+ write(6,*) ' file index = ', file_index
+!
+ open(8,file=calculation(1:l_calc)//file_index//'.g_t',status='unknown',action='write')
+ open(9,file=calculation(1:l_calc)//file_index//'.his',status='unknown',action='write')
+!
+ DO j = 1, 1000
+ read(3,'(a48)',end=305) g_t
+ write(8,'(a48)') g_t
+ END DO
+305 close(8)
+!
+310 read(4,'(a32)',end=390) his
+ IF (his(2:2) .EQ. '-') goto 390
+ write(9,'(a32)') his
+ GOTO 310
+!
+390 close(9)
+ GOTO 300
+!
+! finish
+!
+900 close(1)
+ close(3)
+ close(4)
+ Write(6,*) ' end of *.inp file '
+ Write(6,*) ' '
+ GOTO 200
+ END
\ No newline at end of file
diff --git a/index/index.dsp b/index/index.dsp
new file mode 100755
index 0000000..5e71a76
--- /dev/null
+++ b/index/index.dsp
@@ -0,0 +1,105 @@
+# Microsoft Developer Studio Project File - Name="index" - Package Owner=<4>
+# Microsoft Developer Studio Generated Build File, Format Version 6.00
+# ** DO NOT EDIT **
+
+# TARGTYPE "Win32 (x86) Console Application" 0x0103
+
+CFG=index - Win32 Debug
+!MESSAGE This is not a valid makefile. To build this project using NMAKE,
+!MESSAGE use the Export Makefile command and run
+!MESSAGE
+!MESSAGE NMAKE /f "index.mak".
+!MESSAGE
+!MESSAGE You can specify a configuration when running NMAKE
+!MESSAGE by defining the macro CFG on the command line. For example:
+!MESSAGE
+!MESSAGE NMAKE /f "index.mak" CFG="index - Win32 Debug"
+!MESSAGE
+!MESSAGE Possible choices for configuration are:
+!MESSAGE
+!MESSAGE "index - Win32 Release" (based on "Win32 (x86) Console Application")
+!MESSAGE "index - Win32 Debug" (based on "Win32 (x86) Console Application")
+!MESSAGE
+
+# Begin Project
+# PROP AllowPerConfigDependencies 0
+# PROP Scc_ProjName ""
+# PROP Scc_LocalPath ""
+CPP=cl.exe
+F90=df.exe
+RSC=rc.exe
+
+!IF "$(CFG)" == "index - Win32 Release"
+
+# PROP BASE Use_MFC 0
+# PROP BASE Use_Debug_Libraries 0
+# PROP BASE Output_Dir "Release"
+# PROP BASE Intermediate_Dir "Release"
+# PROP BASE Target_Dir ""
+# PROP Use_MFC 0
+# PROP Use_Debug_Libraries 0
+# PROP Output_Dir "Release"
+# PROP Intermediate_Dir "Release"
+# PROP Target_Dir ""
+# ADD BASE F90 /compile_only /nologo /warn:nofileopt
+# ADD F90 /compile_only /nologo /warn:nofileopt
+# ADD BASE CPP /nologo /W3 /GX /O2 /D "WIN32" /D "NDEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /c
+# ADD CPP /nologo /W3 /GX /O2 /D "WIN32" /D "NDEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /c
+# ADD BASE RSC /l 0x409 /d "NDEBUG"
+# ADD RSC /l 0x409 /d "NDEBUG"
+BSC32=bscmake.exe
+# ADD BASE BSC32 /nologo
+# ADD BSC32 /nologo
+LINK32=link.exe
+# ADD BASE LINK32 kernel32.lib /nologo /subsystem:console /machine:I386
+# ADD LINK32 kernel32.lib /nologo /subsystem:console /machine:I386
+
+!ELSEIF "$(CFG)" == "index - Win32 Debug"
+
+# PROP BASE Use_MFC 0
+# PROP BASE Use_Debug_Libraries 1
+# PROP BASE Output_Dir "Debug"
+# PROP BASE Intermediate_Dir "Debug"
+# PROP BASE Target_Dir ""
+# PROP Use_MFC 0
+# PROP Use_Debug_Libraries 1
+# PROP Output_Dir "Debug"
+# PROP Intermediate_Dir "Debug"
+# PROP Target_Dir ""
+# ADD BASE F90 /check:bounds /compile_only /debug:full /nologo /traceback /warn:argument_checking /warn:nofileopt
+# ADD F90 /check:bounds /compile_only /debug:full /nologo /traceback /warn:argument_checking /warn:nofileopt
+# ADD BASE CPP /nologo /W3 /Gm /GX /ZI /Od /D "WIN32" /D "_DEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /GZ /c
+# ADD CPP /nologo /W3 /Gm /GX /ZI /Od /D "WIN32" /D "_DEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /GZ /c
+# ADD BASE RSC /l 0x409 /d "_DEBUG"
+# ADD RSC /l 0x409 /d "_DEBUG"
+BSC32=bscmake.exe
+# ADD BASE BSC32 /nologo
+# ADD BSC32 /nologo
+LINK32=link.exe
+# ADD BASE LINK32 kernel32.lib /nologo /subsystem:console /debug /machine:I386 /pdbtype:sept
+# ADD LINK32 kernel32.lib /nologo /subsystem:console /incremental:no /debug /machine:I386 /pdbtype:sept
+
+!ENDIF
+
+# Begin Target
+
+# Name "index - Win32 Release"
+# Name "index - Win32 Debug"
+# Begin Group "Source Files"
+
+# PROP Default_Filter "cpp;c;cxx;rc;def;r;odl;idl;hpj;bat;f90;for;f;fpp"
+# Begin Source File
+
+SOURCE=".\index old files.f90"
+# End Source File
+# End Group
+# Begin Group "Header Files"
+
+# PROP Default_Filter "h;hpp;hxx;hm;inl;fi;fd"
+# End Group
+# Begin Group "Resource Files"
+
+# PROP Default_Filter "ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe"
+# End Group
+# End Target
+# End Project
diff --git a/index/index.plg b/index/index.plg
new file mode 100755
index 0000000..7a4cfd3
--- /dev/null
+++ b/index/index.plg
@@ -0,0 +1,19 @@
+
+
+
+Build Log
+
+--------------------Configuration: index - Win32 Release--------------------
+
+Command Lines
+Creating command line "link.exe kernel32.lib /nologo /subsystem:console /incremental:no /pdb:"Release/index.pdb" /machine:I386 /out:"Release/index.exe" ".\Release\index old files.obj" "
+Output Window
+Linking...
+
+
+
+Results
+index.exe - 0 error(s), 0 warning(s)
+
+
+
diff --git a/make spinglass/Debug/DF60.PDB b/make spinglass/Debug/DF60.PDB
new file mode 100755
index 0000000..246f81a
Binary files /dev/null and b/make spinglass/Debug/DF60.PDB differ
diff --git a/make spinglass/Debug/Zero Temperature.obj b/make spinglass/Debug/Zero Temperature.obj
new file mode 100755
index 0000000..e214dbd
Binary files /dev/null and b/make spinglass/Debug/Zero Temperature.obj differ
diff --git a/make spinglass/Debug/make spinglass.pdb b/make spinglass/Debug/make spinglass.pdb
new file mode 100755
index 0000000..f9e4648
Binary files /dev/null and b/make spinglass/Debug/make spinglass.pdb differ
diff --git a/make spinglass/Release/Shortcut to make spinglass.exe.lnk b/make spinglass/Release/Shortcut to make spinglass.exe.lnk
new file mode 100755
index 0000000..bb8c75a
Binary files /dev/null and b/make spinglass/Release/Shortcut to make spinglass.exe.lnk differ
diff --git a/make spinglass/Release/Zero Temperature.obj b/make spinglass/Release/Zero Temperature.obj
new file mode 100755
index 0000000..a4633bc
Binary files /dev/null and b/make spinglass/Release/Zero Temperature.obj differ
diff --git a/make spinglass/Release/field_simulation.obj b/make spinglass/Release/field_simulation.obj
new file mode 100755
index 0000000..5ce8fda
Binary files /dev/null and b/make spinglass/Release/field_simulation.obj differ
diff --git a/make spinglass/Release/make spinglass.exe b/make spinglass/Release/make spinglass.exe
new file mode 100755
index 0000000..a14ac6f
Binary files /dev/null and b/make spinglass/Release/make spinglass.exe differ
diff --git a/make spinglass/Zero Temperature.f90 b/make spinglass/Zero Temperature.f90
new file mode 100755
index 0000000..3a04852
--- /dev/null
+++ b/make spinglass/Zero Temperature.f90
@@ -0,0 +1,385 @@
+! Program to simulate thin film spinglasses.
+!
+! Ge Nieuwenhuys, June, 2002, Written as Ising Metropolis program
+! October 2005, Rewritten as Heisenberg Zero-Temperature
+! October 17, 2005 Bug in direct access file "removes" by
+! oversizing the recordlength
+!
+! October 12: periodic boundary conditions in y- z plane
+! October 14: random number start randomly (based on clock) for
+! batch calculations.
+! October 14: output-file-names are automatically indexed.
+!
+! Spins are located on a fcc lattice
+!
+! nspin number of spins
+! nsp number of spins asked
+! d thickness
+! a lattice constant
+! ah half of lattice constant
+!
+ Use DFPORT ! library only needed for obtaining CPU-time
+ Use DFLIB
+!
+ parameter ( max_spins = 100000 )
+!
+ structure /spin/
+ integer*4 x,y,z
+ real*8 c(3)
+ end structure
+!
+ structure /inter/
+ integer si
+ real*8 val
+ end structure
+!
+ integer*4 j,k,l,m,n, nsp, nspin, nat, steps_per_spin, mh, ix, iy, iz
+ integer*4 l_calc, ifile, iseed
+ record /spin/ s(max_spins)
+ record /inter/ i(100,max_spins)
+ real*8 d, concentration, c, dd(max_spins), rkky, norm, p(3), k_F
+ real*8 b_ext(3), b_ext_K(3), moment, T_glass, mag(3)
+ character*1 answer
+ character*80 calculation, line
+ logical open_inp
+ character*8 dddd, tttt, zone
+ character*4 file_index
+ integer dt(8)
+ real*4 runtime(2), start_time, end_time
+!
+! initialization
+!
+ val = 0
+ a = .407 ! nm, lattice parameter for Au
+ k_F = 12.0 ! 1/nm, Fermi wavevector for Au
+ moment = 2.2 ! mu_B, moment of the impurity spins
+ b_ext(1) = 0.01
+ b_ext(2) = 0.0
+ b_ext(3) = 0.0 ! extrenal field of 100 gauss in x-direction
+ T_glass = 15.0 ! glass_temperature in Kelvin
+ ah = a / 2.0
+ iseed = 1234567
+!
+ Write(6,*) '-----------------------------------------------------'
+ Write(6,*) '| SPIN-GLASS GROUNDSTATE SIMULATION |'
+ Write(6,*) '| Version October 17, 2005 |'
+ Write(6,*) '| |'
+ Write(6,*) '| This program simulates a spin-glass groundstate |'
+ Write(6,*) '| using the method as described by |'
+ Write(6,*) '| Walstedt and Walker, Phys. Rev. B22 (1980) 3816 |'
+ Write(6,*) '| |'
+ Write(6,*) '| The program can be used in batch mode by |'
+ Write(6,*) '| supplying the name of the calculation |'
+ Write(6,*) '| (Enter will start the online mode) |'
+ Write(6,*) '| |'
+ Write(6,*) '| In de batch-mode the parameters will be read from |'
+ Write(6,*) '| .inp |'
+ Write(6,*) '| A comment on the first line and then on |'
+ Write(6,*) '| each line of this file |'
+ Write(6,*) '| number of spins, |'
+ Write(6,*) '| thickness of sample [nm], |'
+ Write(6,*) '| lattice parameter [nm], |'
+ Write(6,*) '| Fermi wave number [nm^-1] |'
+ Write(6,*) '| concentration of spins [at.%], |'
+ Write(6,*) '| glass temperature [K], |'
+ Write(6,*) '| magnetic moment [mu_B], |'
+ Write(6,*) '| external magnetic field (3 components) [tesla] |'
+ Write(6,*) '| number of iterations |'
+ Write(6,*) '| |'
+ Write(6,*) '| can be entered as a |'
+ Write(6,*) '| commandline parameter |'
+ Write(6,*) '-----------------------------------------------------'
+!
+! Ask name of calculation (and read if in batch mode )
+!
+!
+! files :
+!
+ open(9,file='u:\simulations\counter.his',status='old')
+ read(9,*) ifile ! initialize outputfile counter
+!
+ IF ( iargc() .GT. 0 ) THEN
+ call getarg(1, calculation)
+ Write(6,*) ' Calculation taken from commandline > ',calculation
+ ELSE
+777 write(6,778)
+778 format(' '/' Give name of the calculation > ', \)
+ read(5,'(a80)') calculation
+ END IF
+!
+ IF ( calculation(1:1) .NE. ' ' ) THEN
+ l_calc = index( calculation, ' ') - 1
+ open(1,file=calculation(1:l_calc)//'.inp', status='old', action='read', err=777)
+ open(2,file=calculation(1:l_calc)//'.out', status='unknown', action='write')
+ END IF
+!
+ inquire(1,opened=open_inp)
+ IF (open_inp) read(1,'(a80)') comment
+!
+ open(9,file='u:\simulations\counter.sgl',status='old')
+ read(9,*) ifile
+!
+888 IF (open_inp) THEN ! new calculation
+889 ifile = ifile + 1 ! increase outputfile number
+ rewind(9)
+ write(9,*) ifile ! store for next program
+ write(file_index,'(''_'',i3)') ifile ! generate file_name
+ DO j = 2, 4
+ IF (file_index(j:j) .EQ. ' ' ) file_index(j:j) = '0'
+ END DO
+!
+! Although the length of the record /spin/ is 3*4 + 3*8 = 36, the length had to be
+! set to 40, otherwise the writing went wrong (s.c(3) always, except for the last
+! equal to zero
+!
+ open(3,file=calculation(1:l_calc)//file_index//'.sgl',status='unknown', &
+ & access='direct',form='binary',recl=40)
+!
+890 read(1,'(a80)',end=999) line
+ IF ( line(1:1) .EQ. '!' ) THEN
+ write(6,'(a)') line
+ GOTO 890
+ ELSE
+ read(line,*,err=998,end=999) nsp, d, a, k_F, concentration, &
+ & T_glass, moment, b_ext, &
+ & steps_per_spin
+ END IF
+!
+! Initialize randomnumber generator "randomly"
+!
+ call date_and_time( dddd, tttt, zone, dt )
+ DO j = 1, dt(8) ! number milliseconds on the clock
+ dummy = rand(iseed)
+ END DO
+!
+ ELSE
+!
+! Ask size of the system
+!
+1 write(6,2)
+2 format( ' How many spins ? '\)
+ read(5,*,err=1) nsp
+ IF ( nsp .LE. 0 ) STOP ' Programm terminated by operator '
+ IF ( nsp .GT. max_spins ) GOTO 1
+!
+3 write(6,4)
+4 format( ' What thickness [nm] ? '\)
+ read(5,*,err=3) d
+!
+5 write(6,6)
+6 format( ' Which concentration [at.%] ? '\)
+ read(5,*,err=5) concentration
+
+7 write(6,8)
+8 format(' Give the T_glass and the magnetic moment > '\)
+ read(5,*,err=7) T_glass, moment
+ END IF
+!
+! end of getting all parameters
+!
+ c = concentration / 100
+!
+! Calculate the magnetic field energy in Kelvin
+!
+ b_ext_K = moment * (0.927 / 1.38) * b_ext
+!
+ start_time = dtime(runtime) ! record the starttime
+!
+! Calculate the 'rounded' number of spins for a lattice n*m*m for
+! the given concentration
+!
+ n = floor( d / ah )
+ nat = floor( nsp / c )
+ m = floor( sqrt( 2.0 * float(nat) / float(n) ) )
+ mh = m / 2
+ nat = m * m * n / 2
+ nspin = floor( nat * c )
+!
+ write(6,*) ' n = ', n,' m = ', m
+ write(6,*) ' nat = ', nat
+ write(6,*) ' nspin = ', nspin
+!
+ write(2,*) ' n = ', n,' m = ', m
+ write(2,*) ' nat = ', nat
+ write(2,*) ' nspin = ', nspin
+!
+! Place the spins on the lattice
+!
+ nspin = 0
+ DO j = 0, n-1
+ DO k = 0, m-1
+ DO l = 0, m-1
+ IF ( mod(j+k+l,2) .EQ. 0 ) THEN
+ IF ( ran(iseed) .LT. c ) THEN
+ nspin = nspin + 1
+ s(nspin).x = j
+ s(nspin).y = k
+ s(nspin).z = l
+ END IF
+ END IF
+ END DO
+ END DO
+ END DO
+!
+! Calculate the 100 shortest distances
+!
+ write(6,*) ' '
+ write(6,*) ' '
+!
+ DO j = 1, nspin
+ write(6,9) j, char(13)
+9 format(' Considering spin ', i5, a1,\)
+ DO k = 1, nspin
+ IF ( k .EQ. j ) THEN
+ dd(k) = 1e10
+ ELSE
+ ix = s(j).x - s(k).x
+ iy = s(j).y - s(k).y
+ IF ( iy .LT. -mh ) iy = iy + m ! periodic boundary along y-axis
+ IF ( iy .GT. mh ) iy = iy - m
+ iz = s(j).z - s(k).z
+ IF ( iz .LT. -mh ) iz = iz + m ! periodic boundary along z-axis
+ IF ( iz .GT. mh ) iz = iz - m
+ dd(k) = (ix*ix + iy*iy + iz*iz)
+ END IF
+ END DO
+!
+ DO ii = 1, 100
+ dd_min = 1E+10
+ DO k = 1, nspin
+ IF ( dd(k) .LT. dd_min ) THEN
+ l = k
+ dd_min = dd(k)
+ END IF
+ END DO
+ i(ii,j).si = l
+ i(ii,j).val = ah * sqrt(dd(l))
+ dd(l) = 1e10
+ END DO
+!
+! translate distance into interaction strength
+!
+ DO ii = 1, 100
+ i(ii,j).val = rkky( i(ii,j).val, k_F, T_glass )
+ END DO
+!
+ END DO
+!
+ end_time = dtime(runtime)
+ write(6,*) ' Finally, nspin = ', nspin,' in ', end_time - start_time,' seconds '
+ write(2,*) ' Finally, nspin = ', nspin,' in ', end_time - start_time,' seconds '
+!
+! Initialize the spins
+!
+ start_time = dtime(runtime)
+!
+ DO j = 1, nspin
+ DO k = 1, 3
+ s(j).c(k) = 2.0*ran(iseed) - 1.0
+ END DO
+ norm = sqrt( sum( s(j).c * s(j).c ) )
+ s(j).c = s(j).c / norm
+ END DO
+!
+! Calculated the energy and magnetization
+!
+97 e = 0.0
+ mag = 0.0
+!
+ DO j = 1, nspin
+ p = 0.0
+ DO ii = 1, 100
+ p = p + i(ii,j).val * s(i(ii,j).si).c
+ END DO
+ e = e + sum( p * s(j).c ) + sum( b_ext_K * s(j).c )
+ mag = mag + s(j).c
+ END DO
+!
+ end_time = dtime(runtime)
+!
+ Write(6,971) mag / float(nspin)
+ Write(2,971) mag / float(nspin)
+971 format(' Magnetization = ',3F8.4)
+ Write(6,972) e / float(nspin), end_time - start_time
+ Write(2,972) e / float(nspin), end_time - start_time
+972 format( ' Energy = ', E14.4, ' after ', f8.2, ' seconds ' )
+!
+! Now start the serious running
+!
+91 IF ( .NOT. open_inp ) THEN
+98 write(6,99)
+99 format(' Give the number of steps per spin [0: new glass] > '\)
+ read(5,*,err=98) steps_per_spin
+ IF ( steps_per_spin .LE. 0 ) GOTO 777
+ END IF
+!
+ start_time = dtime(runtime) ! record the starttime
+ write(6,*) ' '
+ Write(6,*) ' ' ! to make space for the hashes
+!
+! Now comes the real hard work !!!!!!!!!!!
+!
+ DO mon = 1, steps_per_spin
+ DO j = 1, nspin
+ e = 0.0
+ p = 0.0
+ DO ii = 1, 100
+ p = p + i(ii,j).val * s(i(ii,j).si).c
+ END DO
+ p = p + b_ext_K
+ norm = sqrt( sum( p * p ) )
+ p = p / norm
+ s(j).c = p
+ END DO
+ IF ( mod(mon,100) .EQ. 0 ) idummy = putc('#')
+ END DO
+!
+876 write(6,*) ' '
+ write(6,*) ' '
+ write(3,rec=1) n,m,nspin,a,moment,T_glass
+ write(3,rec=2) concentration,b_ext,steps_per_spin
+ DO ispin = 1, nspin
+ write(3,rec=ispin+2) s(ispin)
+ END DO
+ close(3)
+!
+ IF ( open_inp ) THEN
+ DO j = 1, nspin
+ p = 0.0
+ DO ii = 1, 100
+ p = p + i(ii,j).val * s(i(ii,j).si).c
+ END DO
+ e = e + sum( p * s(j).c ) + sum( b_ext_K * s(j).c )
+ mag = mag + s(j).c
+ END DO
+!
+ end_time = dtime(runtime)
+!
+ Write(6,971) mag / float(nspin)
+ Write(2,971) mag / float(nspin)
+ Write(6,972) e / float(nspin), end_time - start_time
+ Write(2,972) e / float(nspin), end_time - start_time
+ GOTO 888
+ ELSE
+ GOTO 97
+ END IF
+!
+998 STOP 'ERROR IN INPUT FILE '
+!
+999 STOP 'stopped at end of input'
+!
+ END
+!
+!
+ REAL*8 function RKKY(x, k_F, T_glass)
+ real*8 x, k_F, T_glass, xx
+!
+! calculates the RKKY interaction,
+! Te factor of one thousand makes takes care that the RKKY still only
+! has to be multiplied by the glass-temperature.
+!
+ xx = 2.0 * k_F * x
+ rkky = 1.0D+03 * T_glass * ( xx * cos(xx) - sin(xx) ) /(xx*xx*xx*xx)
+!
+ RETURN
+ END
diff --git a/make spinglass/make spinglass.dsp b/make spinglass/make spinglass.dsp
new file mode 100755
index 0000000..1aaa87b
--- /dev/null
+++ b/make spinglass/make spinglass.dsp
@@ -0,0 +1,108 @@
+# Microsoft Developer Studio Project File - Name="make spinglass" - Package Owner=<4>
+# Microsoft Developer Studio Generated Build File, Format Version 6.00
+# ** DO NOT EDIT **
+
+# TARGTYPE "Win32 (x86) Console Application" 0x0103
+
+CFG=make spinglass - Win32 Debug
+!MESSAGE This is not a valid makefile. To build this project using NMAKE,
+!MESSAGE use the Export Makefile command and run
+!MESSAGE
+!MESSAGE NMAKE /f "make spinglass.mak".
+!MESSAGE
+!MESSAGE You can specify a configuration when running NMAKE
+!MESSAGE by defining the macro CFG on the command line. For example:
+!MESSAGE
+!MESSAGE NMAKE /f "make spinglass.mak" CFG="make spinglass - Win32 Debug"
+!MESSAGE
+!MESSAGE Possible choices for configuration are:
+!MESSAGE
+!MESSAGE "make spinglass - Win32 Release" (based on "Win32 (x86) Console Application")
+!MESSAGE "make spinglass - Win32 Debug" (based on "Win32 (x86) Console Application")
+!MESSAGE
+
+# Begin Project
+# PROP AllowPerConfigDependencies 0
+# PROP Scc_ProjName ""
+# PROP Scc_LocalPath ""
+CPP=cl.exe
+F90=df.exe
+RSC=rc.exe
+
+!IF "$(CFG)" == "make spinglass - Win32 Release"
+
+# PROP BASE Use_MFC 0
+# PROP BASE Use_Debug_Libraries 0
+# PROP BASE Output_Dir "Release"
+# PROP BASE Intermediate_Dir "Release"
+# PROP BASE Target_Dir ""
+# PROP Use_MFC 0
+# PROP Use_Debug_Libraries 0
+# PROP Output_Dir "Release"
+# PROP Intermediate_Dir "Release"
+# PROP Ignore_Export_Lib 0
+# PROP Target_Dir ""
+# ADD BASE F90 /compile_only /nologo /warn:nofileopt
+# ADD F90 /compile_only /nologo /warn:nofileopt
+# ADD BASE CPP /nologo /W3 /GX /O2 /D "WIN32" /D "NDEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /c
+# ADD CPP /nologo /W3 /GX /O2 /D "WIN32" /D "NDEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /c
+# ADD BASE RSC /l 0x409 /d "NDEBUG"
+# ADD RSC /l 0x409 /d "NDEBUG"
+BSC32=bscmake.exe
+# ADD BASE BSC32 /nologo
+# ADD BSC32 /nologo
+LINK32=link.exe
+# ADD BASE LINK32 kernel32.lib /nologo /subsystem:console /machine:I386
+# ADD LINK32 kernel32.lib /nologo /subsystem:console /machine:I386 /out:"Release/make_spinglass.exe"
+# SUBTRACT LINK32 /incremental:yes
+
+!ELSEIF "$(CFG)" == "make spinglass - Win32 Debug"
+
+# PROP BASE Use_MFC 0
+# PROP BASE Use_Debug_Libraries 1
+# PROP BASE Output_Dir "Debug"
+# PROP BASE Intermediate_Dir "Debug"
+# PROP BASE Target_Dir ""
+# PROP Use_MFC 0
+# PROP Use_Debug_Libraries 1
+# PROP Output_Dir "Debug"
+# PROP Intermediate_Dir "Debug"
+# PROP Ignore_Export_Lib 0
+# PROP Target_Dir ""
+# ADD BASE F90 /check:bounds /compile_only /debug:full /nologo /traceback /warn:argument_checking /warn:nofileopt
+# ADD F90 /check:bounds /compile_only /debug:full /nologo /traceback /warn:argument_checking /warn:nofileopt
+# ADD BASE CPP /nologo /W3 /Gm /GX /ZI /Od /D "WIN32" /D "_DEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /GZ /c
+# ADD CPP /nologo /W3 /Gm /GX /ZI /Od /D "WIN32" /D "_DEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /GZ /c
+# ADD BASE RSC /l 0x409 /d "_DEBUG"
+# ADD RSC /l 0x409 /d "_DEBUG"
+BSC32=bscmake.exe
+# ADD BASE BSC32 /nologo
+# ADD BSC32 /nologo
+LINK32=link.exe
+# ADD BASE LINK32 kernel32.lib /nologo /subsystem:console /debug /machine:I386 /pdbtype:sept
+# ADD LINK32 kernel32.lib /nologo /subsystem:console /incremental:no /debug /machine:I386 /out:"Debug/make_spinglass.exe" /pdbtype:sept
+
+!ENDIF
+
+# Begin Target
+
+# Name "make spinglass - Win32 Release"
+# Name "make spinglass - Win32 Debug"
+# Begin Group "Source Files"
+
+# PROP Default_Filter "cpp;c;cxx;rc;def;r;odl;idl;hpj;bat;f90;for;f;fpp"
+# Begin Source File
+
+SOURCE=".\Zero Temperature.f90"
+# End Source File
+# End Group
+# Begin Group "Header Files"
+
+# PROP Default_Filter "h;hpp;hxx;hm;inl;fi;fd"
+# End Group
+# Begin Group "Resource Files"
+
+# PROP Default_Filter "ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe"
+# End Group
+# End Target
+# End Project
diff --git a/make spinglass/make spinglass.plg b/make spinglass/make spinglass.plg
new file mode 100755
index 0000000..8acc600
--- /dev/null
+++ b/make spinglass/make spinglass.plg
@@ -0,0 +1,19 @@
+
+
+
+Build Log
+
+--------------------Configuration: make spinglass - Win32 Debug--------------------
+
+Command Lines
+Creating command line "link.exe kernel32.lib /nologo /subsystem:console /incremental:no /pdb:"Debug/make_spinglass.pdb" /debug /machine:I386 /out:"Debug/make_spinglass.exe" /pdbtype:sept ".\Debug\Zero Temperature.obj" "
+Output Window
+Linking...
+
+
+
+Results
+make_spinglass.exe - 0 error(s), 0 warning(s)
+
+
+
diff --git a/monte_carlo.dsw b/monte_carlo.dsw
new file mode 100755
index 0000000..adf0761
--- /dev/null
+++ b/monte_carlo.dsw
@@ -0,0 +1,89 @@
+Microsoft Developer Studio Workspace File, Format Version 6.00
+# WARNING: DO NOT EDIT OR DELETE THIS WORKSPACE FILE!
+
+###############################################################################
+
+Project: "dipole field calculation"=".\dipole field calculation\dipole field calculation.dsp" - Package Owner=<4>
+
+Package=<5>
+{{{
+}}}
+
+Package=<4>
+{{{
+}}}
+
+###############################################################################
+
+Project: "dynamics"=.\dynamics\dynamics.dsp - Package Owner=<4>
+
+Package=<5>
+{{{
+}}}
+
+Package=<4>
+{{{
+}}}
+
+###############################################################################
+
+Project: "index"=.\index\index.dsp - Package Owner=<4>
+
+Package=<5>
+{{{
+}}}
+
+Package=<4>
+{{{
+}}}
+
+###############################################################################
+
+Project: "make spinglass"=".\make spinglass\make spinglass.dsp" - Package Owner=<4>
+
+Package=<5>
+{{{
+}}}
+
+Package=<4>
+{{{
+}}}
+
+###############################################################################
+
+Project: "test"=.\test\test.dsp - Package Owner=<4>
+
+Package=<5>
+{{{
+}}}
+
+Package=<4>
+{{{
+}}}
+
+###############################################################################
+
+Project: "to_plot"=.\to_plot\to_plot.dsp - Package Owner=<4>
+
+Package=<5>
+{{{
+}}}
+
+Package=<4>
+{{{
+}}}
+
+###############################################################################
+
+Global:
+
+Package=<5>
+{{{
+}}}
+
+Package=<3>
+{{{
+}}}
+
+###############################################################################
+
diff --git a/monte_carlo.ncb b/monte_carlo.ncb
new file mode 100755
index 0000000..d586eda
Binary files /dev/null and b/monte_carlo.ncb differ
diff --git a/monte_carlo.opt b/monte_carlo.opt
new file mode 100755
index 0000000..5b4f4a0
Binary files /dev/null and b/monte_carlo.opt differ
diff --git a/test/Debug/DF60.PDB b/test/Debug/DF60.PDB
new file mode 100755
index 0000000..4ce103c
Binary files /dev/null and b/test/Debug/DF60.PDB differ
diff --git a/test/test-1.f90 b/test/test-1.f90
new file mode 100755
index 0000000..5cedad4
--- /dev/null
+++ b/test/test-1.f90
@@ -0,0 +1,91 @@
+!----------------------------------------------------------------------------------------
+! test for precession, ge nieuwenhuys, Oktober 25, 2005
+!----------------------------------------------------------------------------------------
+!
+! The initial muon direction is given by the unit vector emu, the magnetic field
+! is b (vector) and the observer (positroncounter) is supposed to be in the
+! direction of the unit-vector observer.
+!
+ USE DFLIB
+!
+ real*8 b(3), eb(3), emu(3), gyro(3), azi(3), length_b, length_gyro, length_azi
+ real*8 observer(3), amplitude, help(3), cos_phase, phase
+ real*8 length, saclar_pruduct, length_vector_product
+
+1 write(6,2)
+2 format(' Give initial muon > '\)
+ read(5,*,err=1) emu
+ emu = emu / length(emu)
+3 write(6,4)
+4 format(' Give field > '\)
+ read(5,*,err=3) b
+ length_b = length( b )
+ IF ( length_b .LT. 0.001 ) GOTO 1
+5 write(6,6)
+6 format(' Give observer > '\)
+ read(5,*,err=5) observer
+ IF ( length(observer) .LT. 0.01 ) GOTO 3
+ observer = observer / length( observer)
+!
+ eb = b / length_b ! unit vector along b
+ length_azi = scalar_product( emu, eb ) ! length of azimuth of cone
+ azi = length_azi * eb ! vector with length parallel to field
+ length_gyro = length_vector_product( emu, eb ) ! length of the radius of base circle
+ gyro = emu - azi ! that radius as vector at t = 0
+ amplitude = length_vector_product( observer, eb ) ! the projection of the radius on the
+! direction of the observer.
+!
+! Spin will rotate around the field, and start -of course- from the initial spin direction.
+! So the radial vector at t=0 equals emu - azi, while the amplitude is the length
+! of the vector product observer, field. That means that the cosine of the phase angle
+! is the projection of gyro on the observer direction divided by the amplitude.
+!
+ IF ( abs( amplitude ) .GT. 1D-09 ) THEN
+ cos_phase = scalar_product( gyro, observer ) / amplitude
+ ELSE
+ cos_phase = 1.0
+ END IF
+!
+ phase = acos( cos_phase )
+!
+ write(6,'(a,3f10.5)') ' emu = ', emu
+ write(6,'(a,3f10.5)') ' eb = ', eb
+ write(6,'(a,1f10.5)') ' length_azi = ', length_azi
+ write(6,'(a,3f10.5)') ' azi = ', azi
+ write(6,'(a,1f10.5)') ' length_gyro= ', length_gyro
+ write(6,'(a,3f10.5)') ' gyro = ', gyro
+ write(6,'(a,1f10.5)') ' amplitude = ', amplitude
+ write(6,'(a,1f10.5)') ' cos_phase = ', cos_phase
+ write(6,'(a,1f10.5)') ' phase = ', phase
+ write(6,*) ' '
+!
+ GOTO 5
+!
+ END
+!
+ real*8 FUNCTION length( v )
+ real*8 v(3)
+ length = sqrt( sum( v * v ) )
+ RETURN
+ END
+!
+ real*8 FUNCTION scalar_product( v, w )
+ real*8 v(3), w(3)
+ scalar_product = sum( v * w )
+ RETURN
+ END
+!
+ real*8 FUNCTION length_vector_product( v, w )
+ real*8 v(3), w(3), vp(3), length
+ call vector_product( vp, v, w )
+ length_vector_product = length( vp )
+ RETURN
+ END
+!
+ SUBROUTINE vector_product( vp, v, w )
+ real*8 v(3), w(3), vp(3)
+ vp(1) = v(2) * w(3) - v(3) * w(2)
+ vp(2) = v(3) * w(1) - v(1) * w(3)
+ vp(3) = v(1) * w(2) - v(2) * w(1)
+ RETURN
+ END
\ No newline at end of file
diff --git a/test/test-2.f90 b/test/test-2.f90
new file mode 100755
index 0000000..7516f35
--- /dev/null
+++ b/test/test-2.f90
@@ -0,0 +1,62 @@
+!----------------------------------------------------------------------------------------
+! test for precession, ge nieuwenhuys, Oktober 25, 2005
+!----------------------------------------------------------------------------------------
+!
+! The initial muon direction is given by the unit vector emu, the magnetic field
+! is b (vector) and the observer (positroncounter) is supposed to be in the
+! direction of the unit-vector observer.
+!
+ USE DFLIB
+!
+ real*8 b(3), eb(3), emu(3), gyro(3), azi(3), length_b, length_gyro, length_azi
+ real*8 observer(3), amplitude, help(3), cos_phase, phase, g_t_x, g_t_y, g_t_z, omega
+ real*8 length, saclar_pruduct, length_vector_product
+
+1 write(6,2)
+2 format(' Give initial muon > '\)
+ read(5,*,err=1) emu
+ emu = emu / length(emu)
+3 write(6,4)
+4 format(' Give field > '\)
+ read(5,*,err=3) b
+ length_b = length( b )
+ IF ( length_b .LT. 0.001 ) GOTO 1
+!
+ open(1,file='e:\simulations\rotatie.dat',status='unknown')
+ eb = b / length_b
+ omega = 135. * 6.28 * length_b
+!
+ DO i = 1, 100
+!
+ t = float(i-1) * 0.01
+!
+ c = cos( omega * t )
+ s = sin( omega * t )
+!
+ g_t_x = ( c + eb(1)*eb(1)*(1-c) ) * emu(1) + &
+ & ( eb(1)*eb(2)*(1-c)-eb(3)*s ) * emu(2) + &
+ & ( eb(2)*s+eb(1)*eb(3)*(1-c) ) * emu(3)
+!
+ g_t_y = ( eb(3)*s+eb(1)*eb(2)*(1-c) ) * emu(1) + &
+ & ( c + eb(2)*eb(2) * (1-c) ) * emu(2) + &
+ & ( -eb(1)*s+eb(2)*eb(3)*(1-c)) * emu(3)
+!
+ g_t_z = ( -eb(2)*s+eb(1)*eb(3)*(1-c)) * emu(1) + &
+ & ( eb(1)*s+eb(2)*eb(3)*(1-c)) * emu(2) + &
+ & ( c+eb(3)*eb(3)*(1-c) ) * emu(3)
+ !
+ write(1,'(4E18.6)') t, g_t_x, g_t_y, g_t_z
+!
+ END DO
+!
+ Close(1)
+ GOTO 3
+!
+ END
+!
+ real*8 FUNCTION length( v )
+ real*8 v(3)
+ length = sqrt( sum( v * v ) )
+ RETURN
+ END
+!
diff --git a/test/test.dsp b/test/test.dsp
new file mode 100755
index 0000000..a060c29
--- /dev/null
+++ b/test/test.dsp
@@ -0,0 +1,118 @@
+# Microsoft Developer Studio Project File - Name="test" - Package Owner=<4>
+# Microsoft Developer Studio Generated Build File, Format Version 6.00
+# ** DO NOT EDIT **
+
+# TARGTYPE "Win32 (x86) Console Application" 0x0103
+
+CFG=test - Win32 Debug
+!MESSAGE This is not a valid makefile. To build this project using NMAKE,
+!MESSAGE use the Export Makefile command and run
+!MESSAGE
+!MESSAGE NMAKE /f "test.mak".
+!MESSAGE
+!MESSAGE You can specify a configuration when running NMAKE
+!MESSAGE by defining the macro CFG on the command line. For example:
+!MESSAGE
+!MESSAGE NMAKE /f "test.mak" CFG="test - Win32 Debug"
+!MESSAGE
+!MESSAGE Possible choices for configuration are:
+!MESSAGE
+!MESSAGE "test - Win32 Release" (based on "Win32 (x86) Console Application")
+!MESSAGE "test - Win32 Debug" (based on "Win32 (x86) Console Application")
+!MESSAGE
+
+# Begin Project
+# PROP AllowPerConfigDependencies 0
+# PROP Scc_ProjName ""
+# PROP Scc_LocalPath ""
+CPP=cl.exe
+F90=df.exe
+RSC=rc.exe
+
+!IF "$(CFG)" == "test - Win32 Release"
+
+# PROP BASE Use_MFC 0
+# PROP BASE Use_Debug_Libraries 0
+# PROP BASE Output_Dir "Release"
+# PROP BASE Intermediate_Dir "Release"
+# PROP BASE Target_Dir ""
+# PROP Use_MFC 0
+# PROP Use_Debug_Libraries 0
+# PROP Output_Dir "Release"
+# PROP Intermediate_Dir "Release"
+# PROP Target_Dir ""
+# ADD BASE F90 /compile_only /nologo /warn:nofileopt
+# ADD F90 /compile_only /nologo /warn:nofileopt
+# ADD BASE CPP /nologo /W3 /GX /O2 /D "WIN32" /D "NDEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /c
+# ADD CPP /nologo /W3 /GX /O2 /D "WIN32" /D "NDEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /c
+# ADD BASE RSC /l 0x409 /d "NDEBUG"
+# ADD RSC /l 0x409 /d "NDEBUG"
+BSC32=bscmake.exe
+# ADD BASE BSC32 /nologo
+# ADD BSC32 /nologo
+LINK32=link.exe
+# ADD BASE LINK32 kernel32.lib /nologo /subsystem:console /machine:I386
+# ADD LINK32 kernel32.lib /nologo /subsystem:console /machine:I386
+
+!ELSEIF "$(CFG)" == "test - Win32 Debug"
+
+# PROP BASE Use_MFC 0
+# PROP BASE Use_Debug_Libraries 1
+# PROP BASE Output_Dir "Debug"
+# PROP BASE Intermediate_Dir "Debug"
+# PROP BASE Target_Dir ""
+# PROP Use_MFC 0
+# PROP Use_Debug_Libraries 1
+# PROP Output_Dir "Debug"
+# PROP Intermediate_Dir "Debug"
+# PROP Target_Dir ""
+# ADD BASE F90 /check:bounds /compile_only /debug:full /nologo /traceback /warn:argument_checking /warn:nofileopt
+# ADD F90 /check:bounds /compile_only /debug:full /nologo /traceback /warn:argument_checking /warn:nofileopt
+# ADD BASE CPP /nologo /W3 /Gm /GX /ZI /Od /D "WIN32" /D "_DEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /GZ /c
+# ADD CPP /nologo /W3 /Gm /GX /ZI /Od /D "WIN32" /D "_DEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /GZ /c
+# ADD BASE RSC /l 0x409 /d "_DEBUG"
+# ADD RSC /l 0x409 /d "_DEBUG"
+BSC32=bscmake.exe
+# ADD BASE BSC32 /nologo
+# ADD BSC32 /nologo
+LINK32=link.exe
+# ADD BASE LINK32 kernel32.lib /nologo /subsystem:console /debug /machine:I386 /pdbtype:sept
+# ADD LINK32 kernel32.lib /nologo /subsystem:console /incremental:no /debug /machine:I386 /pdbtype:sept
+
+!ENDIF
+
+# Begin Target
+
+# Name "test - Win32 Release"
+# Name "test - Win32 Debug"
+# Begin Group "Source Files"
+
+# PROP Default_Filter "cpp;c;cxx;rc;def;r;odl;idl;hpj;bat;f90;for;f;fpp"
+# Begin Source File
+
+SOURCE=".\test-1.f90"
+
+!IF "$(CFG)" == "test - Win32 Release"
+
+# PROP Exclude_From_Build 1
+
+!ELSEIF "$(CFG)" == "test - Win32 Debug"
+
+!ENDIF
+
+# End Source File
+# Begin Source File
+
+SOURCE=".\test-2.f90"
+# End Source File
+# End Group
+# Begin Group "Header Files"
+
+# PROP Default_Filter "h;hpp;hxx;hm;inl;fi;fd"
+# End Group
+# Begin Group "Resource Files"
+
+# PROP Default_Filter "ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe"
+# End Group
+# End Target
+# End Project
diff --git a/test/test.f90 b/test/test.f90
new file mode 100755
index 0000000..ac0db76
--- /dev/null
+++ b/test/test.f90
@@ -0,0 +1,24 @@
+! test.f90
+!
+! FUNCTIONS:
+! test - Entry point of console application.
+!
+! Example of displaying 'Hello World' at execution time.
+!
+
+!****************************************************************************
+!
+! PROGRAM: test
+!
+! PURPOSE: Entry point for 'Hello World' sample console application.
+!
+!****************************************************************************
+
+ program test
+
+ implicit none
+
+ print *, 'Hello World'
+
+ end program test
+
diff --git a/test/test.plg b/test/test.plg
new file mode 100755
index 0000000..a68bd1a
--- /dev/null
+++ b/test/test.plg
@@ -0,0 +1,19 @@
+
+
+
+Build Log
+
+--------------------Configuration: test - Win32 Release--------------------
+
+Command Lines
+Creating command line "link.exe kernel32.lib /nologo /subsystem:console /incremental:no /pdb:"Release/test.pdb" /machine:I386 /out:"Release/test.exe" ".\Release\test-2.obj" "
+Output Window
+Linking...
+
+
+
+Results
+test.exe - 0 error(s), 0 warning(s)
+
+
+
diff --git a/thinfilm/0.4 2.0 b/thinfilm/0.4 2.0
new file mode 100755
index 0000000..250cbb5
--- /dev/null
+++ b/thinfilm/0.4 2.0
@@ -0,0 +1,60 @@
+
+ ------------------------------------------------------------------------------
+ u:\simulations\output_5.txt
+ ------------------------------------------------------------------------------
+ Calculation started 2005- 9-29 at 15:32
+ ------------------------------------------------------------------------------
+
+ sample = 5.0 nanometer thick, and 30.0 nanometer wide.
+ concentration = 2.0 at. %
+ anisotropy = F
+ number of muons = 5000
+ number of spins = 1
+ average field = -0.597E-06 -0.612E-06 0.236E-06 tesla
+ second moment = 0.584E-05 0.472E-05 0.363E-05 tesla
+ cpu_time = -0.313E-01 seconds
+
+ ------------------------------------------------------------------------------
+ u:\simulations\output_5.txt
+ ------------------------------------------------------------------------------
+ Calculation started 2005- 9-29 at 15:32
+ ------------------------------------------------------------------------------
+
+ sample = 5.0 nanometer thick, and 50.0 nanometer wide.
+ concentration = 2.0 at. %
+ anisotropy = F
+ number of muons = 5000
+ number of spins = 4
+ average field = -0.448E-05 0.175E-07 0.759E-06 tesla
+ second moment = 0.317E-04 0.236E-04 0.221E-04 tesla
+ cpu_time = -0.156E-01 seconds
+
+ ------------------------------------------------------------------------------
+ u:\simulations\output_5.txt
+ ------------------------------------------------------------------------------
+ Calculation started 2005- 9-29 at 15:32
+ ------------------------------------------------------------------------------
+
+ sample = 5.0 nanometer thick, and 15.0 nanometer wide.
+ concentration = 2.0 at. %
+ anisotropy = F
+ number of muons = 5000
+ number of spins = 3
+ average field = 0.780E-05 0.826E-06 -0.503E-05 tesla
+ second moment = 0.228E-04 0.128E-04 0.193E-04 tesla
+ cpu_time = 0.000E+00 seconds
+
+ ------------------------------------------------------------------------------
+ u:\simulations\output_5.txt
+ ------------------------------------------------------------------------------
+ Calculation started 2005- 9-29 at 15:32
+ ------------------------------------------------------------------------------
+
+ sample = 5.0 nanometer thick, and 10.0 nanometer wide.
+ concentration = 2.0 at. %
+ anisotropy = F
+ number of muons = 5000
+ number of spins = 1
+ average field = -0.323E-05 -0.223E-05 0.446E-05 tesla
+ second moment = 0.483E-05 0.404E-05 0.583E-05 tesla
+ cpu_time = 0.000E+00 seconds
diff --git a/thinfilm/field_calculation.f90 b/thinfilm/field_calculation.f90
new file mode 100755
index 0000000..4f26133
--- /dev/null
+++ b/thinfilm/field_calculation.f90
@@ -0,0 +1,644 @@
+! Program to calculate dipolar fields in spinglasses,
+! their distribution and the depolarization of the muon
+!
+! Ge Nieuwenhuys, March, September, October 2005
+!
+! October 12: periodic boundary conditions in y- z plane
+! October 14: random number start randomly (based on clock) for
+! batch calculations.
+! October 14: output-file-names are automatically indexed.
+! October 17: oversized the recordlength of the direct-accessfile for
+! unknown, but apparently essential reasons.
+!
+! Spins are located on a fcc lattice
+!
+! nspin number of spins
+! nsp number of spins asked
+! d thickness
+! a lattice constant
+! ah half of lattice constant
+!
+ Use DFPORT ! library only needed for obtaining CPU-time
+ Use DFLIB
+!
+! Structure to store the position (as lattice site-indexes)
+! and the direction-cosines of each spin.
+!
+ structure /spin/
+ integer*4 x,y,z
+ real*8 dir(3)
+ end structure
+!
+! Declarations, maximumnumber of spins: max_spins, maxd is the maximum number of
+! unitcell-distance for which the spin in included in the calculation
+!
+ parameter( max_spins = 3000000, & ! maximum number of magnetic moments
+ & gyro = 135.5, & ! gyromagnetic ratio of muon
+ & twpi = 6.2831, & ! two times Pi
+ & radius = 2.0, & ! maxinum distance [nm] for
+ ! the dipole-field will be calculated
+ & range = 10.0, & ! maximum absolute value of the field expected
+ & mrange = 4000, & ! range of the integer histograms
+ & nrange = 80 ) ! range of the normalized histograms
+!
+ character*10 dddd, tttt, zone
+ character*4 file_index
+ integer*4 dt(8), ifile, l_calc
+ character*80 comment, calculation, line
+ logical in_open, out_open, g_t_open, his_open, sgl, sgl_open
+ integer*4 j,k,l,m,n, nsp, nspin, nat, id, ihist(3,-mrange:mrange)
+ integer*4 iseed, maxfield, minfield, ihis, ibin, nd1, nd2, kd, ld, mh
+ record /spin/ s(max_spins)
+ real*8 d, concentration, c, dd(max_spins), w, depth1, depth2
+ real*8 px(max_spins),py(max_spins), pz(max_spins)
+ real*8 b(3), factor, moment, help, r_3, r_5, r(3), p_r, sq_3, h(3)
+ real*8 fraction, norm, aver_b(3), sigma_b(3), delta(3), anisotropy, b_ext(3)
+ real*8 g_t(3,0:999), omega, b_abs, b_sq, ca_sq, his, radiussq
+ real*4 runtime(2), start_time, end_time
+ real*8 eb(3), emu(3), cc, ss, theta, phi
+!
+ Write(6,*) ' '
+ Write(6,*) ' ---------------------------------------------------------------------'
+ Write(6,*) ' | Program field-calculation of muons due to random static spins |'
+ Write(6,*) ' | Version of October 31, 2005 |'
+ Write(6,*) ' | |'
+ Write(6,*) ' | Input can also be read from an input file that should be named |'
+ Write(6,*) ' | .inp and contain: |'
+ Write(6,*) ' | |'
+ Write(6,*) ' | ext. field(3) ,thickness, width, c, number_of_muons, |'
+ Write(6,*) ' | lattice-constant [nm], magnetic moment [mu_B], |'
+ Write(6,*) ' | initial-muon-direction(theta, phi)[degrees], |'
+ Write(6,*) ' | (muon-positions from) depth1, (to) depth2 [nm], |'
+ Write(6,*) ' | anisotropy [isotropic=1, planar <1, axial >1 |'
+ Write(6,*) ' | (neg: ferromagnetic along the |'
+ Write(6,*) ' | x - axis (anisotropy = -1.0) |'
+ Write(6,*) ' | y - axis (anisotropy = -2.0) |'
+ Write(6,*) ' | z - axis (anisotropy = -3.0) |'
+ Write(6,*) ' | |'
+ Write(6,*) ' | O R |'
+ Write(6,*) ' | |'
+ Write(6,*) ' | name of the .sgl file produced by |'
+ Write(6,*) ' | MAKE SPINGLASS (starting on the first position), |'
+ Write(6,*) ' | number_of_muons, |'
+ Write(6,*) ' | initial-muon-direction(theta, phi)[degrees], |'
+ Write(6,*) ' | (muon-positions from) depth1, (to) depth2 [nm], |'
+ Write(6,*) ' | |'
+ write(6,*) ' | Lines starting with ! (first position) are treated as comments. |'
+ Write(6,*) ' | can be issued as a commandline parameter |'
+ Write(6,*) ' ---------------------------------------------------------------------'
+!
+! files :
+!
+ open(9,file='\simulations\counter.his',status='old')
+ read(9,*) ifile ! initialize outputfile counter
+!
+! write(6,*) ' iargc = ', iargc()
+ IF ( iargc() .GT. 0 ) THEN
+ call getarg(1, calculation)
+ Write(6,*) ' Calculation taken from commandline > ',calculation
+ ELSE
+200 write(6,201)
+201 format(' '/' Give name of the calculation > ', \)
+ read(5,'(a60)') calculation
+ END IF
+!
+ l_calc = index( calculation, ' ') - 1
+!
+ IF ( l_calc .GT. 0 ) THEN
+ open(1,file=calculation(1:l_calc)//'.inp',status='old',action='read',err=200 )
+ open(2,file=calculation(1:l_calc)//'.out',status='unknown',action='write')
+!
+
+ END IF
+!
+ inquire(1, opened = in_open )
+ inquire(2, opened = out_open )
+!
+! initialization of randomumber generator
+!
+ iseed = 1234567
+!
+! Get eventually other values from the iput file
+!
+111 IF (in_open) THEN
+!
+! Read everything from the input file, one line per calculation
+!
+ ifile = ifile + 1 ! increase outputfile number
+ rewind(9)
+ write(9,*) ifile ! store for next program
+ write(file_index,'(''_'',i3)') ifile ! generate file_name
+ DO j = 2, 4
+ IF (file_index(j:j) .EQ. ' ' ) file_index(j:j) = '0'
+ END DO
+!
+
+ open(3,file=calculation(1:l_calc)//file_index//'.g_t',status='unknown',action='write')
+ open(4,file=calculation(1:l_calc)//file_index//'.his',status='unknown',action='write')
+!
+ inquire(3, opened = g_t_open )
+ inquire(4, opened = his_open )
+!
+112 read(1,'(a80)',end=999) line
+ IF ( ( line(1:1) .GE. 'a' .AND. line(1:1) .LE. 'z' ) .OR. &
+ & ( line(1:1) .GE. 'A' .AND. line(1:1) .LE. 'Z' ) ) THEN
+ l = index( line, ' ') - 1
+ write(6,*) line(1:l)
+ open(7,file=line(1:l)//'.sgl',status='old', &
+ & access='direct',form='binary',recl=40,action='read',err=998)
+ read(line(l+1:80),*,err=998,end=999) n_site, theta, phi, depth1, depth2
+ ELSE
+ IF ( line(1:1) .EQ. '!' ) THEN
+ write(2,'(a)') line
+ GOTO 112
+ ELSE
+ read(line,*,err=998,end=999) a, moment, b_ext, d, w, concentration, &
+ & n_site, theta, phi, depth1, depth2, anisotropy
+ END IF
+ END IF
+!
+! Initialize randomnumber generator "randomly"
+!
+ call date_and_time( dddd, tttt, zone, dt )
+ DO i = 1, dt(8) ! number milliseconds on the clock
+ dummy = rand(iseed)
+ END DO
+!
+ ELSE
+!
+! put standard values in the case of on-line calculation
+! for the lattice (4 nm), moment (2 uB), external field (0,0,0) and
+! initial_muon_spin in y-direction
+!
+!
+ a = 0.4 ! Assume 0.4 nanometer
+ moment = 2.0 ! Assume 2 Bohrmagneton per spin
+ b_ext = 0.0 ! No external field
+ emu = 0.0
+ emu(2) = 1.0 ! initial muon direction along y-axis
+!
+!
+! Ask size of the system
+!
+3 write(6,4)
+4 format( ' What thickness [nm] (0=stop) ? '\)
+ read(5,*,err=3) d
+ IF ( d .LT. 0.0 ) GOTO 3
+ IF ( d .EQ. 0.0 ) THEN
+ Write(6,*) ' '
+ STOP ' program terminated by operator'
+ END IF
+!
+5 write(6,6)
+6 format( ' What width [nm] ? '\)
+ read(5,*,err=5) w
+ IF ( w .LE. 0.0 ) GOTO 5
+ depth1 = 0.0
+ depth2 = w
+!
+7 write(6,8)
+8 format( ' Which concentration [at.%] ? '\)
+ read(5,*,err=7) concentration
+ IF ( concentration .LE. 0.0 ) GOTO 7
+!
+! Ask for the anisotropy.
+! The random value of the direction cosin in the x-direction is multiplied
+! by anisotropy before normalization
+!
+9 write(6,10)
+10 format( ' The random value of the direction cosin in the x-direction'/ &
+& ' is multiplied by anisotropy before normalization'/ &
+& ' Anisotropy [isotrope == 1] ? '\)
+ read(5,*,err=9) anisotropy
+!
+20 write(6,21)
+21 format( ' Give value of the external field (x=perp to film,'/ &
+& ' y=initial_muon > '\)
+ read(5,*,err=20) b_ext
+!
+ END IF ! end reading from input file / keyboard
+!
+!----------------------------------------------------------------------------------------
+! Start calculation
+!----------------------------------------------------------------------------------------
+ call date_and_time( dddd, tttt, zone, dt )
+!
+! If a spinglass has been simulated by MAKE SPINGLASS, then
+! the .sgl file will be read, ELSE a random
+! glass will be generated here.
+!
+ inquire(7, opened = sgl_open )
+!
+ IF ( sgl_open ) THEN ! spin glass has been made
+ read(7,rec=1) n,m,nspin,a,moment,T_glass
+ read(7,rec=2) concentration,b_ext,steps_per_spin
+ DO ispin = 1, nspin
+ read(7,rec=ispin+2) s(ispin)
+ END DO
+ close(7)
+!
+ ELSE ! spin glass has NOT been made
+!
+ c = concentration / 100.0
+!
+! Calculate the 'rounded' number of spins for a lattice n*m*m for
+! the given concentration.
+! n is the number of atoms (half unitcells) perpendicular
+! to the layer (== x-direction).
+! m is the size of the layer in the y- ad z-direction
+!
+ n = floor(2.0 * d / a ) + 2
+ m = floor(2.0 * w / a ) + 2
+ nat = m * m * n / 2
+ nspin = floor( nat * c )
+!
+ IF (nspin .GE. max_spins ) THEN
+ Write(6,*) ' '
+ Write(6,*) ' Too many spins: ', nspin
+ IF ( out_open ) Write(2,*) ' Too many spins: ', nspin
+ GOTO 111
+ END IF
+!
+! Place the spins randomly on the fcc-lattice
+! Run over a whole simple cubic lattice in steps
+! of half of the fcc-unitcell.
+! Then take care of the fcc-structure and
+! decide whether or not to place a spin.
+!
+ nspin = 0
+!
+ DO j = 0, n-1
+ DO k = 0, m-1
+ DO l = 0, m-1
+ IF ( mod(j+k+l,2) .EQ. 0 ) THEN ! This takes care of the fcc structure.
+ IF ( ran(iseed) .LT. c ) THEN
+ nspin = nspin + 1
+ s(nspin).x = j
+ s(nspin).y = k
+ s(nspin).z = l
+
+ IF (anisotropy .GE. 0.0 ) THEN
+!
+! Give the spin an arbitrary direction
+!
+ DO i = 1, 3
+ h(i) = 2.0D+00 * ran(iseed) - 1.0D+00
+ END DO
+!
+! The anisotropy is taken care off by
+! multiplying the direction cosine in
+! the x-direction with ANOSOTROPY
+! before normalizing the direction cosines.
+!
+ h(1) = anisotropy * h(1)
+ norm = sum( h * h )
+ h = h / sqrt( norm )
+ ELSE
+ h = 0.0
+ h(-int(anisotropy)) = 1.0
+ END IF
+
+ s(nspin).dir = h
+!
+ END IF
+ END IF
+ END DO
+ END DO
+ END DO
+!
+! The sample has been grown now.
+!
+ Write(6,*) ' '
+ Write(6,*) 'The sample has been grown, calculation can start'
+ Write(6,*) ' '
+!
+ END IF ! Of reading ,calculation>.sgl or
+! growing magnetic structure
+!
+! Now start the serious work.
+!
+! Use half of the lattice parameter as unit of length
+!
+ ah = a / 2.0
+!
+! help for periodic boundary conditions
+!
+ mh = m / 2
+!
+! the maximum distance squared in units of ah:
+!
+ radiussq = radius * radius / ( ah * ah )
+!
+! Calculate factor to translate to the correct dimensions.
+!
+! factor is ( mu_o / 4 Pi ) * moment * mu_B / ( ah^3 )
+! -- ALL in MKS units --
+! so that the "field" can be calculated as
+! 1/r^5 ( 3 * (s.dir *** r) * r - r^2 s.dir ),
+! where s.dir is the unit vector to the direction of the magnetic moment,
+! and *** stands for the dot-product.
+!
+ factor = 1D-07 * moment * 9.2740019D-24 / ( ah*ah*ah * 1D-27 )
+!
+! see where the muons should go
+!
+ nd1 = floor( depth1 / ah )
+ nd2 = floor( depth2 / ah )
+ IF ( mod( nd1 , 2 ) .EQ. 0 ) nd1 = nd1 + 1 ! nd1 should be odd
+ IF ( nd2 .LT. nd1 + 1 ) nd2 = nd1 + 1
+!
+! calculate unit vector along the initial muon-spin direction
+!
+ emu(1) = sin( twpi * theta / 360.0 ) * cos( twpi * phi / 360.0)
+ emu(2) = sin( twpi * theta / 360.0 ) * sin( twpi * phi / 360.0)
+ emu(3) = cos( twpi * theta / 360.0 )
+!
+! Ask the number of sites to calculated, about 10,000 is reasonable
+!
+ IF ( .NOT. in_open ) THEN ! read keyboard if no input file
+!
+ write(6,*) ' total number of muon-sites :', (m-1)*(m-1)*(nd2-nd1+1) / 8
+ write(6,*) ' '
+11 write(6,12)
+12 format(' Give number of sites to be calculated > ' $)
+ read(5,*,err=11) n_site
+!
+ END IF ! of reading keyboard
+!
+ fraction = dble( float(n_site) / float( (m-1)*(m-1)*(nd2-nd1+1)/8 ))
+!
+! make some space
+!
+ Write(6,*) ' '
+ Write(6,*) ' '
+!
+ start_time = dtime(runtime) ! record the starttime
+!
+! Initialize the averages
+!
+ ib = 0 ! index of field calculation
+ aver_b = 0 ! average of the field
+ sigma_b = 0 ! average of the field squared
+ hist = 0 ! histograms
+ g_t = 0.0 ! initialize the line
+!
+! Assume the muon to be in the center of the fcc-cube
+!
+ DO j = nd1, nd2, 2
+ DO k = 1, m-1, 2
+ DO l = 1, m-1, 2
+!
+! These do-loops run over all sites, which is probably too much (time consuming)
+! Therefore select randomly sufficient (see above) fraction of
+! the possible muon sites and calculate the dipolar field.
+!
+ IF ( ran(iseed) .LT. fraction ) THEN
+!
+! Calculate the field by running over all spins.
+! In calculating the mutual distance, periodic boundaryconditions are applied
+! in the y- and z-direction, but NOT in the x-direction, since that is supposed
+! perpendicular to the film
+!
+! The field is only calculated when the distance is smaller then radius
+!
+ b = 0
+!
+ DO i = 1, nspin
+ r(1) = dble(float(j-s(i).x))
+ kd = k - s(i).y
+ IF ( kd .LT. -mh ) kd = kd + m ! periodic boundary condition
+ IF ( kd .GT. mh ) kd = kd - m ! periodic boundary condition
+ r(2) = dble(float(kd))
+ ld = l - s(i).z
+ IF ( ld .LT. -mh ) ld = ld + m ! periodic boundary condition
+ IF ( ld .GT. mh ) ld = ld - m ! periodic boundary condition
+ r(3) = dble(float(ld))
+ r_2 = sum( r * r )
+!
+ IF ( r_2 .LE. radiussq ) THEN ! skip calculation if distance is too large
+ help = sqrt( r_2 )
+ r_3 = r_2 * help
+ r_5 = r_2 * r_3
+ h = s(i).dir
+ p_r = sum( h * r )
+ b = b + ( 3.0D+00 * p_r * r - r_2 * h ) / r_5
+ END IF
+!
+ END DO
+!
+ ib = ib + 1 ! count the sites calculated.
+ b = factor * b ! get correct dimensions
+ aver_b = aver_b + b ! add the field to the averages
+ sigma_b = sigma_b + b*b
+!
+!
+! Count for histograms
+!
+ DO ih = 1, 3
+ ival = int( float(mrange) * b(ih) / range + 0.5D+00 )
+ IF ( abs(ival) .LE. mrange ) ihist(ih,ival) = ihist(ih,ival) + 1
+ END DO
+!
+ b = b + b_ext ! add external field
+ b_sq = sum( b * b ) ! square of the field
+ b_abs = sqrt( b_sq ) ! absolute value
+ eb = b / b_abs ! unit vector
+ omega = gyro * twpi * b_abs ! precession frequency
+!
+! Calculate the rotation of the muonspin for 1000 time-steps.
+! The contribution to the asymmetry equals the components of the temporal
+! muonspin, assuming the counters to be forward-backward, left-right ,and up-down,
+! respectively.
+!
+ DO it = 0, 999
+ t = 1.0D-02 * dble(float(it))
+ cc = cos( omega * t )
+ ss = sin( omega * t )
+!
+ g_t(1,it) = g_t(1,it) + &
+ & ( cc+eb(1)*eb(1)*(1-cc)) * emu(1) + &
+ & ( -eb(3)*ss+eb(1)*eb(2)*(1-cc)) * emu(2) + &
+ & ( eb(2)*ss+eb(1)*eb(3)*(1-cc)) * emu(3)
+!
+ g_t(2,it) = g_t(2,it) + &
+ & ( eb(3)*ss+eb(1)*eb(2)*(1-cc)) * emu(1) + &
+ & ( cc+eb(2)*eb(2)*(1-cc)) * emu(2) + &
+ & ( -eb(1)*ss+eb(2)*eb(3)*(1-cc)) * emu(3)
+!
+ g_t(3,it) = g_t(3,it) + &
+ & ( -eb(2)*ss+eb(1)*eb(3)*(1-cc)) * emu(1) + &
+ & ( eb(1)*ss+eb(2)*eb(3)*(1-cc)) * emu(2) + &
+ & ( cc+eb(3)*eb(3)*(1-cc)) * emu(3)
+!
+ END DO
+!
+ IF ( mod(ib,1000) .EQ. 0 ) idummy = putc('#')
+!
+ END IF ! decision on fraction of muon sites
+ END DO
+ END DO
+ END DO ! l, k, j loops
+!
+! Average over all calculaled sites.
+!
+ norm = dble( float(ib))
+ aver_b = aver_b / norm
+ sigma_b = sqrt( (sigma_b - aver_b * aver_b ) / norm )
+ delta = gyro * sigma_b
+ g_t = g_t / norm
+!
+! Renormalize histograms
+!
+ IF ( his_open ) THEN ! Should the histogram be calculated ??
+ Write(4,*) '-------------------------------------------------------'
+!
+! Check whether the maximum calculated field exceeds the range
+!
+ IF ( ihist(1,-mrange) .EQ. 0 .AND. ihist(1,mrange) .EQ. 0 .AND. &
+ & ihist(2,-mrange) .EQ. 0 .AND. ihist(2,mrange) .EQ. 0 .AND. &
+ & ihist(3,-mrange) .EQ. 0 .AND. ihist(3,mrange) .EQ. 0 ) THEN
+!
+! determine the range of fields found
+!
+ DO j = 1, 3
+ DO k = -mrange, mrange
+ IF ( ihist(j, k) .GT. 0 ) maxfield = k
+ IF ( ihist(j,-k) .GT. 0 ) minfield = -k
+ END DO
+!
+! adjust binning of histogram and write values
+!
+ ibin = (maxfield - minfield) / nrange + 1
+ x = float(minfield) * range / float(mrange)
+ step = range * float(ibin) / float(mrange)
+!
+ write(6,*) ' The field histogram vaues are: '
+ write(6,*) minfield, maxfield, ibin, x, step
+!
+ DO i = minfield, maxfield, ibin
+ ihis = 0
+ DO k = 0, ibin-1
+ ihis = ihis + ihist(j,i+k)
+ END DO
+ his = float(ihis) / norm
+ Write(4,'(2E16.6)') x, his
+ x = x + step
+ END DO
+!
+ Write(4,*) ' '
+ END DO
+!
+ ELSE
+ Write(4,*) ' Fields exceed the maximum field for histogram calculation '
+ END IF
+ END IF ! Histogram calculation
+!
+ end_time = dtime(runtime)
+!
+ write(6,*) ' '
+ write(2,100) comment(1:73),(dt(j),j=1,3),(dt(j),j=5,8)
+ write(6,101) n*ah, m*ah
+ write(6,301) nd1*ah, nd2*ah
+ write(6,102) concentration
+ write(6,103) anisotropy, int(-anisotropy)
+ write(6,104) n_site
+ write(6,304) theta, phi
+ write(6,105) nspin
+ write(6,106) aver_b
+ write(6,107) sigma_b
+ write(6,108) delta
+ write(6,308) b_ext
+ write(6,109) end_time - start_time
+!
+! Look whether data have to be written to file
+!
+ IF ( out_open ) THEN
+ write(2,100) comment(1:73),(dt(j),j=1,3),(dt(j),j=5,8)
+ write(2,101) n*ah, m*ah
+ write(2,301) nd1*ah, nd2*ah
+ write(2,102) concentration
+ write(2,103) anisotropy, int(-anisotropy)
+ write(2,104) n_site
+ write(2,304) theta, phi
+ write(2,105) nspin
+ write(2,106) aver_b
+ write(2,107) sigma_b
+ write(2,108) delta
+ write(2,308) b_ext
+ write(2,109) end_time - start_time
+ END IF
+!
+100 format(' '/' ',73('-')/' ',a73/' ',73('-')/ &
+& ' Calculation started ',i5,'-',i2,'-',i2, &
+& ' at ',2(i2,':'),i2,'.',i3/' ',73('-')/' ')
+101 format(' sample = ', F6.1, ' nanometer thick, and ', F6.1, ' nanometer wide.')
+102 format(' concentration = ', F12.1, ' at. %')
+103 format(' anisotropy = ', E12.3,' (int) ',I2)
+104 format(' number of muons = ', I12)
+105 format(' number of spins = ', I12)
+106 format(' average field = ', 3E12.3,' tesla')
+107 format(' second moment = ', 3E12.3,' tesla')
+108 format(' corres. delta = ', 3E12.3,' 1/microseconde')
+109 format(' cpu_time = ', E12.3, ' seconds')
+308 format(' ext. field = ', 3E12.3,' tesla')
+301 format(' penetration from = ', F6.1,' to ',F6.1' nanometer.')
+304 format(' initial muon spin, theta = ',f6.2,' phi = ', f6.2)
+!
+! Write G(t) if the file is open
+!
+500 IF ( g_t_open ) THEN
+!
+ DO k = 0, 999
+ write(3,'(3E20.6)') (g_t(id,k),id=1,3) ! output
+ END DO
+!
+ END IF
+!
+! Go back to read new parameters
+!
+ GOTO 111
+!
+! On error in input_file
+!
+998 Write(6,*) ' '
+ Write(6,*) ' There is an error in the input file. '
+ IF ( out_open ) Write(2,*) ' There is an error in the input file. '
+!
+999 IF ( in_open ) close(1)
+ IF ( out_open ) close(2)
+ IF ( g_t_open ) close(3)
+ IF ( his_open ) close(4)
+ END
+!
+! End of program
+!-------------------------------------------------------------------------------------------
+!
+! Functions and Subroutines
+!
+!-------------------------------------------------------------------------------------------
+ real*8 FUNCTION length( v )
+ real*8 v(3)
+ length = sqrt( sum( v * v ) )
+ RETURN
+ END
+!
+ real*8 FUNCTION scalar_product( v, w )
+ real*8 v(3), w(3)
+ scalar_product = sum( v * w )
+ RETURN
+ END
+!
+ real*8 FUNCTION length_vector_product( v, w )
+ real*8 v(3), w(3), vp(3), length
+ call vector_product( vp, v, w )
+ length_vector_product = length( vp )
+ RETURN
+ END
+!
+ SUBROUTINE vector_product( vp, v, w )
+ real*8 v(3), w(3), vp(3)
+ vp(1) = v(2) * w(3) - v(3) * w(2)
+ vp(2) = v(3) * w(1) - v(1) * w(3)
+ vp(3) = v(1) * w(2) - v(2) * w(1)
+ RETURN
+ END
diff --git a/thinfilm/field_calculation_GaAs.f90 b/thinfilm/field_calculation_GaAs.f90
new file mode 100755
index 0000000..4da7a6a
--- /dev/null
+++ b/thinfilm/field_calculation_GaAs.f90
@@ -0,0 +1,716 @@
+! Program to calculate dipolar fields in spinglasses,
+! their distribution and the depolarization of the muon
+!
+! Ge Nieuwenhuys, March, September, October 2005
+!
+! October 12: periodic boundary conditions in y- z plane
+! October 14: random number start randomly (based on clock) for
+! batch calculations.
+! October 14: output-file-names are automatically indexed.
+! October 17: oversized the recordlength of the direct-accessfile for
+! unknown, but apparently essential reasons.
+!
+! Spins are located on a fcc lattice
+!
+! nspin number of spins
+! nsp number of spins asked
+! d thickness
+! a lattice constant
+! ah half of lattice constant
+!
+ Use DFPORT ! library only needed for obtaining CPU-time
+ Use DFLIB
+!
+! Structure to store the position (as lattice site-indexes)
+! and the direction-cosines of each spin.
+!
+ structure /spin/
+ integer*4 x,y,z
+ real*8 dir(3)
+ end structure
+!
+! Declarations, maximumnumber of spins: max_spins, maxd is the maximum number of
+! unitcell-distance for which the spin in included in the calculation
+!
+ parameter( max_spins = 3000000, & ! maximum number of magnetic moments
+ & gyro = 135.5, & ! gyromagnetic ratio of muon
+ & twpi = 6.2831, & ! two times Pi
+ & radius = 2.0, & ! maxinum distance [nm] for
+ ! the dipole-field will be calculated
+ & range = 10.0, & ! maximum absolute value of the field expected
+ & mrange = 4000, & ! range of the integer histograms
+ & nrange = 80 ) ! range of the normalized histograms
+!
+ character*10 dddd, tttt, zone
+ character*4 file_index
+ integer*4 dt(8), ifile, l_calc, bond
+ character*80 comment, calculation, line
+ logical in_open, out_open, g_t_open, his_open, sgl, sgl_open
+ integer*4 j,k,l,m,n, nsp, nspin, nat, id, ihist(3,-mrange:mrange)
+ integer*4 iseed, maxfield, minfield, ihis, ibin, nd1, nd2, kd, ld, mh
+ record /spin/ s(max_spins)
+ real*8 d, concentration, c, dd(max_spins), w, depth1, depth2
+ real*8 px(max_spins),py(max_spins), pz(max_spins)
+ real*8 b(3), factor, moment, help, r_3, r_5, r(3), p_r, sq_3, h(3)
+ real*8 fraction, norm, aver_b(3), sigma_b(3), delta(3), anisotropy, b_ext(3)
+ real*8 g_t(3,0:999), omega, b_abs, b_sq, ca_sq, his, radiussq
+ real*4 runtime(2), start_time, end_time
+ real*8 eb(3), emu(3), cc, ss, theta, phi
+!
+ real*8 gaas(16,3)
+!
+! Coordinates of bond sites
+!
+ gaas(1,1) = 1
+ gaas(1,2) = 1
+ gaas(1,3) = 7
+
+ gaas(2,1) = 3
+ gaas(2,2) = 3
+ gaas(2,3) = 7
+
+ gaas(3,1) = 1
+ gaas(3,2) = 3
+ gaas(3,3) = 5
+
+ gaas(4,1) = 3
+ gaas(4,2) = 1
+ gaas(4,3) = 5
+
+ gaas(5,1) = 7
+ gaas(5,2) = 7
+ gaas(5,3) = 7
+
+ gaas(6,1) = 5
+ gaas(6,2) = 5
+ gaas(6,3) = 7
+
+ gaas(7,1) = 7
+ gaas(7,2) = 5
+ gaas(7,3) = 5
+
+ gaas(8,1) = 5
+ gaas(8,2) = 7
+ gaas(8,3) = 5
+
+ gaas(9,1) = 3
+ gaas(9,2) = 7
+ gaas(9,3) = 3
+
+ gaas(10,1) = 1
+ gaas(10,2) = 5
+ gaas(10,3) = 3
+
+ gaas(11,1) = 1
+ gaas(11,2) = 7
+ gaas(11,3) = 1
+
+ gaas(12,1) = 3
+ gaas(12,2) = 5
+ gaas(12,3) = 1
+
+ gaas(13,1) = 7
+ gaas(13,2) = 3
+ gaas(13,3) = 3
+
+ gaas(14,1) = 5
+ gaas(14,2) = 1
+ gaas(14,3) = 3
+
+ gaas(15,1) = 7
+ gaas(15,2) = 1
+ gaas(15,3) = 1
+
+ gaas(16,1) = 5
+ gaas(16,2) = 3
+ gaas(16,3) = 1
+!
+ gaas = gaas / 4.0 ! concert to units of half lattice constant
+!
+ Write(6,*) ' '
+ Write(6,*) ' ---------------------------------------------------------------------'
+ Write(6,*) ' | Program field-calculation of muons due to random static spins |'
+ Write(6,*) ' | Version of October 31, 2005 |'
+ Write(6,*) ' | |'
+ Write(6,*) ' | Input can also be read from an input file that should be named |'
+ Write(6,*) ' | .inp and contain: |'
+ Write(6,*) ' | |'
+ Write(6,*) ' | ext. field(3) ,thickness, width, c, number_of_muons, |'
+ Write(6,*) ' | lattice-constant [nm], magnetic moment [mu_B], |'
+ Write(6,*) ' | initial-muon-direction(theta, phi)[degrees], |'
+ Write(6,*) ' | (muon-positions from) depth1, (to) depth2 [nm], |'
+ Write(6,*) ' | anisotropy [isotropic=1, planar <1, axial >1 |'
+ Write(6,*) ' | (neg: ferromagnetic along the |'
+ Write(6,*) ' | x - axis (anisotropy = -1.0) |'
+ Write(6,*) ' | y - axis (anisotropy = -2.0) |'
+ Write(6,*) ' | z - axis (anisotropy = -3.0) |'
+ Write(6,*) ' | |'
+ Write(6,*) ' | O R |'
+ Write(6,*) ' | |'
+ Write(6,*) ' | name of the .sgl file produced by |'
+ Write(6,*) ' | MAKE SPINGLASS (starting on the first position), |'
+ Write(6,*) ' | number_of_muons, |'
+ Write(6,*) ' | initial-muon-direction(theta, phi)[degrees], |'
+ Write(6,*) ' | (muon-positions from) depth1, (to) depth2 [nm], |'
+ Write(6,*) ' | |'
+ write(6,*) ' | Lines starting with ! (first position) are treated as comments. |'
+ Write(6,*) ' | can be issued as a commandline parameter |'
+ Write(6,*) ' ---------------------------------------------------------------------'
+!
+! files :
+!
+ open(9,file='\simulations\counter.his',status='old')
+ read(9,*) ifile ! initialize outputfile counter
+!
+! write(6,*) ' iargc = ', iargc()
+ IF ( iargc() .GT. 0 ) THEN
+ call getarg(1, calculation)
+ Write(6,*) ' Calculation taken from commandline > ',calculation
+ ELSE
+200 write(6,201)
+201 format(' '/' Give name of the calculation > ', \)
+ read(5,'(a60)') calculation
+ END IF
+!
+ l_calc = index( calculation, ' ') - 1
+!
+ IF ( l_calc .GT. 0 ) THEN
+ open(1,file=calculation(1:l_calc)//'.inp',status='old',action='read',err=200 )
+ open(2,file=calculation(1:l_calc)//'.out',status='unknown',action='write')
+!
+
+ END IF
+!
+ inquire(1, opened = in_open )
+ inquire(2, opened = out_open )
+!
+! initialization of randomumber generator
+!
+ iseed = 1234567
+!
+! Get eventually other values from the iput file
+!
+111 IF (in_open) THEN
+!
+! Read everything from the input file, one line per calculation
+!
+ ifile = ifile + 1 ! increase outputfile number
+ rewind(9)
+ write(9,*) ifile ! store for next program
+ write(file_index,'(''_'',i3)') ifile ! generate file_name
+ DO j = 2, 4
+ IF (file_index(j:j) .EQ. ' ' ) file_index(j:j) = '0'
+ END DO
+!
+
+ open(3,file=calculation(1:l_calc)//file_index//'.g_t',status='unknown',action='write')
+ open(4,file=calculation(1:l_calc)//file_index//'.his',status='unknown',action='write')
+!
+ inquire(3, opened = g_t_open )
+ inquire(4, opened = his_open )
+!
+112 read(1,'(a80)',end=999) line
+ IF ( ( line(1:1) .GE. 'a' .AND. line(1:1) .LE. 'z' ) .OR. &
+ & ( line(1:1) .GE. 'A' .AND. line(1:1) .LE. 'Z' ) ) THEN
+ l = index( line, ' ') - 1
+ write(6,*) line(1:l)
+ open(7,file=line(1:l)//'.sgl',status='old', &
+ & access='direct',form='binary',recl=40,action='read',err=998)
+ read(line(l+1:80),*,err=998,end=999) n_site, theta, phi, depth1, depth2
+ ELSE
+ IF ( line(1:1) .EQ. '!' ) THEN
+ write(2,'(a)') line
+ GOTO 112
+ ELSE
+ read(line,*,err=998,end=999) a, moment, b_ext, d, w, concentration, &
+ & n_site, theta, phi, depth1, depth2, anisotropy
+ END IF
+ END IF
+!
+! Initialize randomnumber generator "randomly"
+!
+ call date_and_time( dddd, tttt, zone, dt )
+ DO i = 1, dt(8) ! number milliseconds on the clock
+ dummy = rand(iseed)
+ END DO
+!
+ ELSE
+!
+! put standard values in the case of on-line calculation
+! for the lattice (4 nm), moment (2 uB), external field (0,0,0) and
+! initial_muon_spin in y-direction
+!
+!
+ a = 0.4 ! Assume 0.4 nanometer
+ moment = 2.0 ! Assume 2 Bohrmagneton per spin
+ b_ext = 0.0 ! No external field
+ emu = 0.0
+ emu(2) = 1.0 ! initial muon direction along y-axis
+!
+!
+! Ask size of the system
+!
+3 write(6,4)
+4 format( ' What thickness [nm] (0=stop) ? '\)
+ read(5,*,err=3) d
+ IF ( d .LT. 0.0 ) GOTO 3
+ IF ( d .EQ. 0.0 ) THEN
+ Write(6,*) ' '
+ STOP ' program terminated by operator'
+ END IF
+!
+5 write(6,6)
+6 format( ' What width [nm] ? '\)
+ read(5,*,err=5) w
+ IF ( w .LE. 0.0 ) GOTO 5
+ depth1 = 0.0
+ depth2 = w
+!
+7 write(6,8)
+8 format( ' Which concentration [at.%] ? '\)
+ read(5,*,err=7) concentration
+ IF ( concentration .LE. 0.0 ) GOTO 7
+!
+! Ask for the anisotropy.
+! The random value of the direction cosin in the x-direction is multiplied
+! by anisotropy before normalization
+!
+9 write(6,10)
+10 format( ' The random value of the direction cosin in the x-direction'/ &
+& ' is multiplied by anisotropy before normalization'/ &
+& ' Anisotropy [isotrope == 1] ? '\)
+ read(5,*,err=9) anisotropy
+!
+20 write(6,21)
+21 format( ' Give value of the external field (x=perp to film,'/ &
+& ' y=initial_muon > '\)
+ read(5,*,err=20) b_ext
+!
+ END IF ! end reading from input file / keyboard
+!
+!----------------------------------------------------------------------------------------
+! Start calculation
+!----------------------------------------------------------------------------------------
+ call date_and_time( dddd, tttt, zone, dt )
+!
+! If a spinglass has been simulated by MAKE SPINGLASS, then
+! the .sgl file will be read, ELSE a random
+! glass will be generated here.
+!
+ inquire(7, opened = sgl_open )
+!
+ IF ( sgl_open ) THEN ! spin glass has been made
+ read(7,rec=1) n,m,nspin,a,moment,T_glass
+ read(7,rec=2) concentration,b_ext,steps_per_spin
+ DO ispin = 1, nspin
+ read(7,rec=ispin+2) s(ispin)
+ END DO
+ close(7)
+!
+ ELSE ! spin glass has NOT been made
+!
+ c = concentration / 100.0
+!
+! Calculate the 'rounded' number of spins for a lattice n*m*m for
+! the given concentration.
+! n is the number of atoms (half unitcells) perpendicular
+! to the layer (== x-direction).
+! m is the size of the layer in the y- ad z-direction
+!
+ n = floor(2.0 * d / a ) + 2
+ m = floor(2.0 * w / a ) + 2
+ nat = m * m * n / 2
+ nspin = floor( nat * c )
+!
+ IF (nspin .GE. max_spins ) THEN
+ Write(6,*) ' '
+ Write(6,*) ' Too many spins: ', nspin
+ IF ( out_open ) Write(2,*) ' Too many spins: ', nspin
+ GOTO 111
+ END IF
+!
+! Place the spins randomly on the fcc-lattice
+! Run over a whole simple cubic lattice in steps
+! of half of the fcc-unitcell.
+! Then take care of the fcc-structure and
+! decide whether or not to place a spin.
+!
+ nspin = 0
+!
+ DO j = 0, n-1
+ DO k = 0, m-1
+ DO l = 0, m-1
+ IF ( mod(j+k+l,2) .EQ. 0 ) THEN ! This takes care of the fcc structure.
+ IF ( ran(iseed) .LT. c ) THEN
+ nspin = nspin + 1
+ s(nspin).x = j
+ s(nspin).y = k
+ s(nspin).z = l
+
+ IF (anisotropy .GE. 0.0 ) THEN
+!
+! Give the spin an arbitrary direction
+!
+ DO i = 1, 3
+ h(i) = 2.0D+00 * ran(iseed) - 1.0D+00
+ END DO
+!
+! The anisotropy is taken care off by
+! multiplying the direction cosine in
+! the x-direction with ANOSOTROPY
+! before normalizing the direction cosines.
+!
+ h(1) = anisotropy * h(1)
+ norm = sum( h * h )
+ h = h / sqrt( norm )
+ ELSE
+ h = 0.0
+ h(-int(anisotropy)) = 1.0
+ END IF
+
+ s(nspin).dir = h
+!
+ END IF
+ END IF
+ END DO
+ END DO
+ END DO
+!
+! The sample has been grown now.
+!
+ Write(6,*) ' '
+ Write(6,*) 'The sample has been grown, calculation can start'
+ Write(6,*) ' '
+!
+ END IF ! Of reading ,calculation>.sgl or
+! growing magnetic structure
+!
+! Now start the serious work.
+!
+! Use half of the lattice parameter as unit of length
+!
+ ah = a / 2.0
+!
+! help for periodic boundary conditions
+!
+ mh = m / 2
+!
+! the maximum distance squared in units of ah:
+!
+ radiussq = radius * radius / ( ah * ah )
+!
+! Calculate factor to translate to the correct dimensions.
+!
+! factor is ( mu_o / 4 Pi ) * moment * mu_B / ( ah^3 )
+! -- ALL in MKS units --
+! so that the "field" can be calculated as
+! 1/r^5 ( 3 * (s.dir *** r) * r - r^2 s.dir ),
+! where s.dir is the unit vector to the direction of the magnetic moment,
+! and *** stands for the dot-product.
+!
+ factor = 1D-07 * moment * 9.2740019D-24 / ( ah*ah*ah * 1D-27 )
+!
+! see where the muons should go
+!
+ nd1 = floor( depth1 / ah )
+ nd2 = floor( depth2 / ah )
+ IF ( mod( nd1 , 2 ) .EQ. 0 ) nd1 = nd1 + 1 ! nd1 should be odd
+ IF ( nd2 .LT. nd1 + 1 ) nd2 = nd1 + 1
+!
+! calculate unit vector along the initial muon-spin direction
+!
+ emu(1) = sin( twpi * theta / 360.0 ) * cos( twpi * phi / 360.0)
+ emu(2) = sin( twpi * theta / 360.0 ) * sin( twpi * phi / 360.0)
+ emu(3) = cos( twpi * theta / 360.0 )
+!
+! Ask the number of sites to calculated, about 10,000 is reasonable
+!
+ IF ( .NOT. in_open ) THEN ! read keyboard if no input file
+!
+ write(6,*) ' total number of muon-sites :', 2*(m-1)*(m-1)*(nd2-nd1+1)
+ write(6,*) ' '
+11 write(6,12)
+12 format(' Give number of sites to be calculated > ' $)
+ read(5,*,err=11) n_site
+!
+ END IF ! of reading keyboard
+!
+ fraction = dble( float(n_site) / float( 2*(m-1)*(m-1)*(nd2-nd1+1)))
+!
+! make some space
+!
+ Write(6,*) ' '
+ Write(6,*) ' '
+!
+ start_time = dtime(runtime) ! record the starttime
+!
+! Initialize the averages
+!
+ ib = 0 ! index of field calculation
+ aver_b = 0 ! average of the field
+ sigma_b = 0 ! average of the field squared
+ hist = 0 ! histograms
+ g_t = 0.0 ! initialize the line
+!
+! Assume the muon to be in the center of the fcc-cube
+!
+ DO j = nd1-1, nd2-1, 2
+ DO k = 0, m-2, 2
+ DO l = 0, m-2, 2
+ DO bond = 1, 16 ! loop over the 16 bonds
+!
+! These do-loops run over all sites, which is probably too much (time consuming)
+! Therefore select randomly sufficient (see above) fraction of
+! the possible muon sites and calculate the dipolar field.
+!
+ IF ( ran(iseed) .LT. fraction ) THEN
+!
+! Calculate the field by running over all spins.
+! In calculating the mutual distance, periodic boundaryconditions are applied
+! in the y- and z-direction, but NOT in the x-direction, since that is supposed
+! perpendicular to the film
+!
+! The field is only calculated when the distance is smaller then radius
+!
+ b = 0
+!
+ DO i = 1, nspin
+ r(1) = dble(float(j-s(i).x)) + gaas(bond,1)
+ kd = k - s(i).y
+ IF ( kd .LT. -mh ) kd = kd + m ! periodic boundary condition
+ IF ( kd .GT. mh ) kd = kd - m ! periodic boundary condition
+ r(2) = dble(float(kd)) + gaas(bond,2)
+ ld = l - s(i).z
+ IF ( ld .LT. -mh ) ld = ld + m ! periodic boundary condition
+ IF ( ld .GT. mh ) ld = ld - m ! periodic boundary condition
+ r(3) = dble(float(ld)) + gaas(bond,3)
+ r_2 = sum( r * r )
+!
+ IF ( r_2 .LE. radiussq ) THEN ! skip calculation if distance is too large
+ help = sqrt( r_2 )
+ r_3 = r_2 * help
+ r_5 = r_2 * r_3
+ h = s(i).dir
+ p_r = sum( h * r )
+ b = b + ( 3.0D+00 * p_r * r - r_2 * h ) / r_5
+ END IF
+!
+ END DO
+!
+ ib = ib + 1 ! count the sites calculated.
+ b = factor * b ! get correct dimensions
+ aver_b = aver_b + b ! add the field to the averages
+ sigma_b = sigma_b + b*b
+!
+!
+! Count for histograms
+!
+ DO ih = 1, 3
+ ival = int( float(mrange) * b(ih) / range + 0.5D+00 )
+ IF ( abs(ival) .LE. mrange ) ihist(ih,ival) = ihist(ih,ival) + 1
+ END DO
+!
+ b = b + b_ext ! add external field
+ b_sq = sum( b * b ) ! square of the field
+ b_abs = sqrt( b_sq ) ! absolute value
+ eb = b / b_abs ! unit vector
+ omega = gyro * twpi * b_abs ! precession frequency
+!
+! Calculate the rotation of the muonspin for 1000 time-steps.
+! The contribution to the asymmetry equals the components of the temporal
+! muonspin, assuming the counters to be forward-backward, left-right ,and up-down,
+! respectively.
+!
+ DO it = 0, 999
+ t = 1.0D-02 * dble(float(it))
+ cc = cos( omega * t )
+ ss = sin( omega * t )
+!
+ g_t(1,it) = g_t(1,it) + &
+ & ( cc+eb(1)*eb(1)*(1-cc)) * emu(1) + &
+ & ( -eb(3)*ss+eb(1)*eb(2)*(1-cc)) * emu(2) + &
+ & ( eb(2)*ss+eb(1)*eb(3)*(1-cc)) * emu(3)
+!
+ g_t(2,it) = g_t(2,it) + &
+ & ( eb(3)*ss+eb(1)*eb(2)*(1-cc)) * emu(1) + &
+ & ( cc+eb(2)*eb(2)*(1-cc)) * emu(2) + &
+ & ( -eb(1)*ss+eb(2)*eb(3)*(1-cc)) * emu(3)
+!
+ g_t(3,it) = g_t(3,it) + &
+ & ( -eb(2)*ss+eb(1)*eb(3)*(1-cc)) * emu(1) + &
+ & ( eb(1)*ss+eb(2)*eb(3)*(1-cc)) * emu(2) + &
+ & ( cc+eb(3)*eb(3)*(1-cc)) * emu(3)
+!
+ END DO
+!
+ IF ( mod(ib,1000) .EQ. 0 ) idummy = putc('#')
+!
+ END IF ! decision on fraction of muon sites
+ END DO ! over bond loop
+ END DO
+ END DO
+ END DO ! l, k, j loops
+!
+! Average over all calculaled sites.
+!
+ norm = dble( float(ib))
+ aver_b = aver_b / norm
+ sigma_b = sqrt( (sigma_b - aver_b * aver_b ) / norm )
+ delta = gyro * sigma_b
+ g_t = g_t / norm
+!
+! Renormalize histograms
+!
+ IF ( his_open ) THEN ! Should the histogram be calculated ??
+ Write(4,*) '-------------------------------------------------------'
+!
+! Check whether the maximum calculated field exceeds the range
+!
+ IF ( ihist(1,-mrange) .EQ. 0 .AND. ihist(1,mrange) .EQ. 0 .AND. &
+ & ihist(2,-mrange) .EQ. 0 .AND. ihist(2,mrange) .EQ. 0 .AND. &
+ & ihist(3,-mrange) .EQ. 0 .AND. ihist(3,mrange) .EQ. 0 ) THEN
+!
+! determine the range of fields found
+!
+ DO j = 1, 3
+ DO k = -mrange, mrange
+ IF ( ihist(j, k) .GT. 0 ) maxfield = k
+ IF ( ihist(j,-k) .GT. 0 ) minfield = -k
+ END DO
+!
+! adjust binning of histogram and write values
+!
+ ibin = (maxfield - minfield) / nrange + 1
+ x = float(minfield) * range / float(mrange)
+ step = range * float(ibin) / float(mrange)
+!
+ write(6,*) ' The field histogram vaues are: '
+ write(6,*) minfield, maxfield, ibin, x, step
+!
+ DO i = minfield, maxfield, ibin
+ ihis = 0
+ DO k = 0, ibin-1
+ ihis = ihis + ihist(j,i+k)
+ END DO
+ his = float(ihis) / norm
+ Write(4,'(2E16.6)') x, his
+ x = x + step
+ END DO
+!
+ Write(4,*) ' '
+ END DO
+!
+ ELSE
+ Write(4,*) ' Fields exceed the maximum field for histogram calculation '
+ END IF
+ END IF ! Histogram calculation
+!
+ end_time = dtime(runtime)
+!
+ write(6,*) ' '
+ write(2,100) comment(1:73),(dt(j),j=1,3),(dt(j),j=5,8)
+ write(6,101) n*ah, m*ah
+ write(6,301) nd1*ah, nd2*ah
+ write(6,102) concentration
+ write(6,103) anisotropy, int(-anisotropy)
+ write(6,104) n_site
+ write(6,304) theta, phi
+ write(6,105) nspin
+ write(6,106) aver_b
+ write(6,107) sigma_b
+ write(6,108) delta
+ write(6,308) b_ext
+ write(6,109) end_time - start_time
+!
+! Look whether data have to be written to file
+!
+ IF ( out_open ) THEN
+ write(2,100) comment(1:73),(dt(j),j=1,3),(dt(j),j=5,8)
+ write(2,101) n*ah, m*ah
+ write(2,301) nd1*ah, nd2*ah
+ write(2,102) concentration
+ write(2,103) anisotropy, int(-anisotropy)
+ write(2,104) n_site
+ write(2,304) theta, phi
+ write(2,105) nspin
+ write(2,106) aver_b
+ write(2,107) sigma_b
+ write(2,108) delta
+ write(2,308) b_ext
+ write(2,109) end_time - start_time
+ END IF
+!
+100 format(' '/' ',73('-')/' ',a73/' ',73('-')/ &
+& ' Calculation started ',i5,'-',i2,'-',i2, &
+& ' at ',2(i2,':'),i2,'.',i3/' ',73('-')/' ')
+101 format(' sample = ', F6.1, ' nanometer thick, and ', F6.1, ' nanometer wide.')
+102 format(' concentration = ', F12.1, ' at. %')
+103 format(' anisotropy = ', E12.3,' (int) ',I2)
+104 format(' number of muons = ', I12)
+105 format(' number of spins = ', I12)
+106 format(' average field = ', 3E12.3,' tesla')
+107 format(' second moment = ', 3E12.3,' tesla')
+108 format(' corres. delta = ', 3E12.3,' 1/microseconde')
+109 format(' cpu_time = ', E12.3, ' seconds')
+308 format(' ext. field = ', 3E12.3,' tesla')
+301 format(' penetration from = ', F6.1,' to ',F6.1' nanometer.')
+304 format(' initial muon spin, theta = ',f6.2,' phi = ', f6.2)
+!
+! Write G(t) if the file is open
+!
+500 IF ( g_t_open ) THEN
+!
+ DO k = 0, 999
+ write(3,'(3E20.6)') (g_t(id,k),id=1,3) ! output
+ END DO
+!
+ END IF
+!
+! Go back to read new parameters
+!
+ GOTO 111
+!
+! On error in input_file
+!
+998 Write(6,*) ' '
+ Write(6,*) ' There is an error in the input file. '
+ IF ( out_open ) Write(2,*) ' There is an error in the input file. '
+!
+999 IF ( in_open ) close(1)
+ IF ( out_open ) close(2)
+ IF ( g_t_open ) close(3)
+ IF ( his_open ) close(4)
+ END
+!
+! End of program
+!-------------------------------------------------------------------------------------------
+!
+! Functions and Subroutines
+!
+!-------------------------------------------------------------------------------------------
+ real*8 FUNCTION length( v )
+ real*8 v(3)
+ length = sqrt( sum( v * v ) )
+ RETURN
+ END
+!
+ real*8 FUNCTION scalar_product( v, w )
+ real*8 v(3), w(3)
+ scalar_product = sum( v * w )
+ RETURN
+ END
+!
+ real*8 FUNCTION length_vector_product( v, w )
+ real*8 v(3), w(3), vp(3), length
+ call vector_product( vp, v, w )
+ length_vector_product = length( vp )
+ RETURN
+ END
+!
+ SUBROUTINE vector_product( vp, v, w )
+ real*8 v(3), w(3), vp(3)
+ vp(1) = v(2) * w(3) - v(3) * w(2)
+ vp(2) = v(3) * w(1) - v(1) * w(3)
+ vp(3) = v(1) * w(2) - v(2) * w(1)
+ RETURN
+ END
diff --git a/thinfilm/test.f90 b/thinfilm/test.f90
new file mode 100755
index 0000000..0ea7b68
--- /dev/null
+++ b/thinfilm/test.f90
@@ -0,0 +1,3 @@
+ Write(6,*) ' het werkt weer '
+ stop
+ end
diff --git a/thinfilm/test.g_t b/thinfilm/test.g_t
new file mode 100755
index 0000000..e69de29
diff --git a/thinfilm/thinfilm.001 b/thinfilm/thinfilm.001
new file mode 100755
index 0000000..a29852e
--- /dev/null
+++ b/thinfilm/thinfilm.001
@@ -0,0 +1,98 @@
+# Microsoft Developer Studio Project File - Name="thinfilm" - Package Owner=<4>
+# Microsoft Developer Studio Generated Build File, Format Version 5.00
+# ** DO NOT EDIT **
+
+# TARGTYPE "Win32 (x86) Console Application" 0x0103
+
+CFG=thinfilm - Win32 Debug
+!MESSAGE This is not a valid makefile. To build this project using NMAKE,
+!MESSAGE use the Export Makefile command and run
+!MESSAGE
+!MESSAGE NMAKE /f "thinfilm.mak".
+!MESSAGE
+!MESSAGE You can specify a configuration when running NMAKE
+!MESSAGE by defining the macro CFG on the command line. For example:
+!MESSAGE
+!MESSAGE NMAKE /f "thinfilm.mak" CFG="thinfilm - Win32 Debug"
+!MESSAGE
+!MESSAGE Possible choices for configuration are:
+!MESSAGE
+!MESSAGE "thinfilm - Win32 Release" (based on\
+ "Win32 (x86) Console Application")
+!MESSAGE "thinfilm - Win32 Debug" (based on "Win32 (x86) Console Application")
+!MESSAGE
+
+# Begin Project
+# PROP Scc_ProjName ""
+# PROP Scc_LocalPath ""
+F90=df.exe
+RSC=rc.exe
+
+!IF "$(CFG)" == "thinfilm - Win32 Release"
+
+# PROP BASE Use_MFC 0
+# PROP BASE Use_Debug_Libraries 0
+# PROP BASE Output_Dir "Release"
+# PROP BASE Intermediate_Dir "Release"
+# PROP BASE Target_Dir ""
+# PROP Use_MFC 0
+# PROP Use_Debug_Libraries 0
+# PROP Output_Dir "Release"
+# PROP Intermediate_Dir "Release"
+# PROP Target_Dir ""
+# ADD BASE F90 /include:"Release/" /compile_only /nologo /warn:nofileopt
+# ADD F90 /include:"Release/" /compile_only /nologo /warn:nofileopt
+# ADD BASE RSC /l 0x413 /d "NDEBUG"
+# ADD RSC /l 0x413 /d "NDEBUG"
+BSC32=bscmake.exe
+# ADD BASE BSC32 /nologo
+# ADD BSC32 /nologo
+LINK32=link.exe
+# ADD BASE LINK32 kernel32.lib /nologo /subsystem:console /machine:I386
+# ADD LINK32 kernel32.lib /nologo /subsystem:console /machine:I386
+
+!ELSEIF "$(CFG)" == "thinfilm - Win32 Debug"
+
+# PROP BASE Use_MFC 0
+# PROP BASE Use_Debug_Libraries 1
+# PROP BASE Output_Dir "Debug"
+# PROP BASE Intermediate_Dir "Debug"
+# PROP BASE Target_Dir ""
+# PROP Use_MFC 0
+# PROP Use_Debug_Libraries 1
+# PROP Output_Dir "Debug"
+# PROP Intermediate_Dir "Debug"
+# PROP Target_Dir ""
+# ADD BASE F90 /include:"Debug/" /compile_only /nologo /debug:full /optimize:0 /warn:nofileopt
+# ADD F90 /include:"Debug/" /compile_only /nologo /debug:full /optimize:0 /warn:nofileopt
+# ADD BASE RSC /l 0x413 /d "_DEBUG"
+# ADD RSC /l 0x413 /d "_DEBUG"
+BSC32=bscmake.exe
+# ADD BASE BSC32 /nologo
+# ADD BSC32 /nologo
+LINK32=link.exe
+# ADD BASE LINK32 kernel32.lib /nologo /subsystem:console /debug /machine:I386 /pdbtype:sept
+# ADD LINK32 kernel32.lib /nologo /subsystem:console /debug /machine:I386 /pdbtype:sept
+
+!ENDIF
+
+# Begin Target
+
+# Name "thinfilm - Win32 Release"
+# Name "thinfilm - Win32 Debug"
+# Begin Source File
+
+SOURCE=.\field_calculation.f90
+# End Source File
+# Begin Source File
+
+SOURCE=.\field_simulation.f90
+# PROP Exclude_From_Build 1
+# End Source File
+# Begin Source File
+
+SOURCE=.\test.f90
+# PROP Exclude_From_Build 1
+# End Source File
+# End Target
+# End Project
diff --git a/thinfilm/thinfilm.dsp b/thinfilm/thinfilm.dsp
new file mode 100755
index 0000000..5af811e
--- /dev/null
+++ b/thinfilm/thinfilm.dsp
@@ -0,0 +1,106 @@
+# Microsoft Developer Studio Project File - Name="thinfilm" - Package Owner=<4>
+# Microsoft Developer Studio Generated Build File, Format Version 6.00
+# ** DO NOT EDIT **
+
+# TARGTYPE "Win32 (x86) Console Application" 0x0103
+
+CFG=thinfilm - Win32 Debug
+!MESSAGE This is not a valid makefile. To build this project using NMAKE,
+!MESSAGE use the Export Makefile command and run
+!MESSAGE
+!MESSAGE NMAKE /f "thinfilm.mak".
+!MESSAGE
+!MESSAGE You can specify a configuration when running NMAKE
+!MESSAGE by defining the macro CFG on the command line. For example:
+!MESSAGE
+!MESSAGE NMAKE /f "thinfilm.mak" CFG="thinfilm - Win32 Debug"
+!MESSAGE
+!MESSAGE Possible choices for configuration are:
+!MESSAGE
+!MESSAGE "thinfilm - Win32 Release" (based on "Win32 (x86) Console Application")
+!MESSAGE "thinfilm - Win32 Debug" (based on "Win32 (x86) Console Application")
+!MESSAGE
+
+# Begin Project
+# PROP AllowPerConfigDependencies 0
+# PROP Scc_ProjName ""
+# PROP Scc_LocalPath ""
+CPP=cl.exe
+F90=df.exe
+RSC=rc.exe
+
+!IF "$(CFG)" == "thinfilm - Win32 Release"
+
+# PROP BASE Use_MFC 0
+# PROP BASE Use_Debug_Libraries 0
+# PROP BASE Output_Dir "Release"
+# PROP BASE Intermediate_Dir "Release"
+# PROP BASE Target_Dir ""
+# PROP Use_MFC 0
+# PROP Use_Debug_Libraries 0
+# PROP Output_Dir "Release"
+# PROP Intermediate_Dir "Release"
+# PROP Target_Dir ""
+# ADD BASE F90 /compile_only /include:"Release/" /nologo /warn:nofileopt
+# ADD F90 /compile_only /include:"Release/" /nologo /warn:nofileopt
+# ADD BASE RSC /l 0x413 /d "NDEBUG"
+# ADD RSC /l 0x413 /d "NDEBUG"
+BSC32=bscmake.exe
+# ADD BASE BSC32 /nologo
+# ADD BSC32 /nologo
+LINK32=link.exe
+# ADD BASE LINK32 kernel32.lib /nologo /subsystem:console /machine:I386
+# ADD LINK32 kernel32.lib /nologo /subsystem:console /machine:I386
+
+!ELSEIF "$(CFG)" == "thinfilm - Win32 Debug"
+
+# PROP BASE Use_MFC 0
+# PROP BASE Use_Debug_Libraries 1
+# PROP BASE Output_Dir "Debug"
+# PROP BASE Intermediate_Dir "Debug"
+# PROP BASE Target_Dir ""
+# PROP Use_MFC 0
+# PROP Use_Debug_Libraries 1
+# PROP Output_Dir "Debug"
+# PROP Intermediate_Dir "Debug"
+# PROP Target_Dir ""
+# ADD BASE F90 /compile_only /debug:full /include:"Debug/" /nologo /warn:nofileopt
+# ADD F90 /compile_only /debug:full /include:"Debug/" /nologo /warn:nofileopt
+# ADD BASE RSC /l 0x413 /d "_DEBUG"
+# ADD RSC /l 0x413 /d "_DEBUG"
+BSC32=bscmake.exe
+# ADD BASE BSC32 /nologo
+# ADD BSC32 /nologo
+LINK32=link.exe
+# ADD BASE LINK32 kernel32.lib /nologo /subsystem:console /debug /machine:I386 /pdbtype:sept
+# ADD LINK32 kernel32.lib /nologo /subsystem:console /debug /machine:I386 /pdbtype:sept
+
+!ENDIF
+
+# Begin Target
+
+# Name "thinfilm - Win32 Release"
+# Name "thinfilm - Win32 Debug"
+# Begin Source File
+
+SOURCE=.\field_calculation.f90
+
+!IF "$(CFG)" == "thinfilm - Win32 Release"
+
+!ELSEIF "$(CFG)" == "thinfilm - Win32 Debug"
+
+!ENDIF
+
+# End Source File
+# Begin Source File
+
+SOURCE=.\field_simulation.f90
+# PROP Exclude_From_Build 1
+# End Source File
+# Begin Source File
+
+SOURCE=.\test.f90
+# PROP Exclude_From_Build 1
+# End Source File
+# End Target
+# End Project
diff --git a/thinfilm/thinfilm.plg b/thinfilm/thinfilm.plg
new file mode 100755
index 0000000..68eb72f
--- /dev/null
+++ b/thinfilm/thinfilm.plg
@@ -0,0 +1,26 @@
+
+
+
+Build Log
+
+--------------------Configuration: thinfilm - Win32 Release--------------------
+
+Command Lines
+Creating temporary file "C:\DOCUME~1\NIEUWE~1\LOCALS~1\Temp\RSP28A.tmp" with contents
+[
+/compile_only /include:"Release/" /nologo /warn:nofileopt /module:"Release/" /object:"Release/"
+"U:\monte_carlo\thinfilm\field_calculation.f90"
+]
+Creating command line "link.exe kernel32.lib /nologo /subsystem:console /incremental:no /pdb:"Release/thinfilm.pdb" /machine:I386 /out:"Release/thinfilm.exe" .\Release\field_calculation.obj "
+Output Window
+Compiling Fortran...
+U:\monte_carlo\thinfilm\field_calculation.f90
+Linking...
+
+
+
+Results
+thinfilm.exe - 0 error(s), 0 warning(s)
+
+
+
diff --git a/to_plot/Debug/DF60.PDB b/to_plot/Debug/DF60.PDB
new file mode 100755
index 0000000..74f9947
Binary files /dev/null and b/to_plot/Debug/DF60.PDB differ
diff --git a/to_plot/Debug/to_plot.exe b/to_plot/Debug/to_plot.exe
new file mode 100755
index 0000000..7757042
Binary files /dev/null and b/to_plot/Debug/to_plot.exe differ
diff --git a/to_plot/Debug/to_plot.obj b/to_plot/Debug/to_plot.obj
new file mode 100755
index 0000000..be88846
Binary files /dev/null and b/to_plot/Debug/to_plot.obj differ
diff --git a/to_plot/Debug/to_plot.pdb b/to_plot/Debug/to_plot.pdb
new file mode 100755
index 0000000..c02d4ad
Binary files /dev/null and b/to_plot/Debug/to_plot.pdb differ
diff --git a/to_plot/Release/dynamics.obj b/to_plot/Release/dynamics.obj
new file mode 100755
index 0000000..c555118
Binary files /dev/null and b/to_plot/Release/dynamics.obj differ
diff --git a/to_plot/Release/to_plot.exe b/to_plot/Release/to_plot.exe
new file mode 100755
index 0000000..7cee739
Binary files /dev/null and b/to_plot/Release/to_plot.exe differ
diff --git a/to_plot/Release/to_plot.obj b/to_plot/Release/to_plot.obj
new file mode 100755
index 0000000..b6deb40
Binary files /dev/null and b/to_plot/Release/to_plot.obj differ
diff --git a/to_plot/to_plot.dsp b/to_plot/to_plot.dsp
new file mode 100755
index 0000000..491ab11
--- /dev/null
+++ b/to_plot/to_plot.dsp
@@ -0,0 +1,105 @@
+# Microsoft Developer Studio Project File - Name="to_plot" - Package Owner=<4>
+# Microsoft Developer Studio Generated Build File, Format Version 6.00
+# ** DO NOT EDIT **
+
+# TARGTYPE "Win32 (x86) Console Application" 0x0103
+
+CFG=to_plot - Win32 Debug
+!MESSAGE This is not a valid makefile. To build this project using NMAKE,
+!MESSAGE use the Export Makefile command and run
+!MESSAGE
+!MESSAGE NMAKE /f "to_plot.mak".
+!MESSAGE
+!MESSAGE You can specify a configuration when running NMAKE
+!MESSAGE by defining the macro CFG on the command line. For example:
+!MESSAGE
+!MESSAGE NMAKE /f "to_plot.mak" CFG="to_plot - Win32 Debug"
+!MESSAGE
+!MESSAGE Possible choices for configuration are:
+!MESSAGE
+!MESSAGE "to_plot - Win32 Release" (based on "Win32 (x86) Console Application")
+!MESSAGE "to_plot - Win32 Debug" (based on "Win32 (x86) Console Application")
+!MESSAGE
+
+# Begin Project
+# PROP AllowPerConfigDependencies 0
+# PROP Scc_ProjName ""
+# PROP Scc_LocalPath ""
+CPP=cl.exe
+F90=df.exe
+RSC=rc.exe
+
+!IF "$(CFG)" == "to_plot - Win32 Release"
+
+# PROP BASE Use_MFC 0
+# PROP BASE Use_Debug_Libraries 0
+# PROP BASE Output_Dir "Release"
+# PROP BASE Intermediate_Dir "Release"
+# PROP BASE Target_Dir ""
+# PROP Use_MFC 0
+# PROP Use_Debug_Libraries 0
+# PROP Output_Dir "Release"
+# PROP Intermediate_Dir "Release"
+# PROP Target_Dir ""
+# ADD BASE F90 /compile_only /nologo /warn:nofileopt
+# ADD F90 /compile_only /nologo /warn:nofileopt
+# ADD BASE CPP /nologo /W3 /GX /O2 /D "WIN32" /D "NDEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /c
+# ADD CPP /nologo /W3 /GX /O2 /D "WIN32" /D "NDEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /c
+# ADD BASE RSC /l 0x409 /d "NDEBUG"
+# ADD RSC /l 0x409 /d "NDEBUG"
+BSC32=bscmake.exe
+# ADD BASE BSC32 /nologo
+# ADD BSC32 /nologo
+LINK32=link.exe
+# ADD BASE LINK32 kernel32.lib /nologo /subsystem:console /machine:I386
+# ADD LINK32 kernel32.lib /nologo /subsystem:console /machine:I386
+
+!ELSEIF "$(CFG)" == "to_plot - Win32 Debug"
+
+# PROP BASE Use_MFC 0
+# PROP BASE Use_Debug_Libraries 1
+# PROP BASE Output_Dir "Debug"
+# PROP BASE Intermediate_Dir "Debug"
+# PROP BASE Target_Dir ""
+# PROP Use_MFC 0
+# PROP Use_Debug_Libraries 1
+# PROP Output_Dir "Debug"
+# PROP Intermediate_Dir "Debug"
+# PROP Target_Dir ""
+# ADD BASE F90 /check:bounds /compile_only /dbglibs /debug:full /nologo /traceback /warn:argument_checking /warn:nofileopt
+# ADD F90 /check:bounds /compile_only /dbglibs /debug:full /nologo /traceback /warn:argument_checking /warn:nofileopt
+# ADD BASE CPP /nologo /W3 /Gm /GX /ZI /Od /D "WIN32" /D "_DEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /GZ /c
+# ADD CPP /nologo /W3 /Gm /GX /ZI /Od /D "WIN32" /D "_DEBUG" /D "_CONSOLE" /D "_MBCS" /YX /FD /GZ /c
+# ADD BASE RSC /l 0x409 /d "_DEBUG"
+# ADD RSC /l 0x409 /d "_DEBUG"
+BSC32=bscmake.exe
+# ADD BASE BSC32 /nologo
+# ADD BSC32 /nologo
+LINK32=link.exe
+# ADD BASE LINK32 kernel32.lib /nologo /subsystem:console /debug /machine:I386 /pdbtype:sept
+# ADD LINK32 kernel32.lib /nologo /subsystem:console /incremental:no /debug /machine:I386 /pdbtype:sept
+
+!ENDIF
+
+# Begin Target
+
+# Name "to_plot - Win32 Release"
+# Name "to_plot - Win32 Debug"
+# Begin Group "Source Files"
+
+# PROP Default_Filter "cpp;c;cxx;rc;def;r;odl;idl;hpj;bat;f90;for;f;fpp"
+# Begin Source File
+
+SOURCE=.\to_plot.f90
+# End Source File
+# End Group
+# Begin Group "Header Files"
+
+# PROP Default_Filter "h;hpp;hxx;hm;inl;fi;fd"
+# End Group
+# Begin Group "Resource Files"
+
+# PROP Default_Filter "ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe"
+# End Group
+# End Target
+# End Project
diff --git a/to_plot/to_plot.f90 b/to_plot/to_plot.f90
new file mode 100755
index 0000000..780b4c6
--- /dev/null
+++ b/to_plot/to_plot.f90
@@ -0,0 +1,88 @@
+! to_plot.f90
+!
+! FUNCTIONS:
+! to_plot - Entry point of console application.
+!
+
+!****************************************************************************
+!
+! PROGRAM: to_plot
+!
+! PURPOSE: To put several *g_t files in a multicolumn file and in the ZF
+! cases correct for the fact that all simulations were done
+! with 50 degree phase angle (should have been 0 for ZF)
+!
+!****************************************************************************
+
+ program to_plot
+
+ implicit none
+
+ ! Variables
+
+ integer*4 i,j,k,n,is,l,max_l
+ real*4 x(40,1000), y(40,1000), xi, yi, yj
+ character*512 out
+ character*80 filename, file_out
+ character*1 y_n
+ logical*4 ZF
+
+ ! Body of to_plot
+
+1 write(6,2)
+2 format(' Give simulation numbers > '$)
+ read(5,*,err=1) is, k, n
+3 write(6,4)
+4 format(' Zero-Field ? '$)
+ read(5,'(a1)') y_n
+ ZF = ( y_n .EQ. 'y' .OR. y_n .EQ. 'Y' )
+5 write(6,6)
+6 format(' Give output file name > '$)
+ read(5,'(a80)') file_out
+
+ max_l = 0
+ x = 0.0
+ y = 0.0
+
+ DO i = k, n
+ write(filename,50) is, i
+50 format('u:\simulations\dynamics-',i2,'_',i3,'.g_t')
+ write(6,*) filename
+ open(1,file=filename,status='old',err=55)
+ write(out( 14*(i-k)+1:14*(i-k+1) ),51) is,i
+51 format(' time ',i2,'_',i3)
+
+ l = 0
+
+ read(1,*) xi, yi, yj
+
+
+ DO WHILE( xi .LT. 5.0 .AND. (.NOT. Eof(1) ) )
+ l = l + 1
+ IF ( ZF ) yi = sqrt( yi*yi + yj*yj )
+ x(i-k+1,l) = xi
+ y(i-k+1,l) = yi
+ read(1,*) xi, yi, yj
+ END DO
+
+ IF ( max_l .LT. l ) max_l = l
+ close(1)
+
+55 END DO
+
+ open(2,file='u:\simulations\'//file_out,status='new')
+ write( 2, '(a)' ) out(1:14*(n-k+1))
+
+ DO j = 1, max_l
+ write(out,'(35(f6.3,f8.3))') ( (x(i,j),y(i,j)),i=1,n-k+1 )
+ DO i = 1, n-k+1
+ IF ( x(i,j) .LT. 0.0 ) out( 14*(i-1)+1:14*i ) = ' '
+ END DO
+ write(2,'(a)') out(1:14*(n-k+1))
+ END DO
+
+ close(2)
+ goto 1
+
+ end program to_plot
+
diff --git a/to_plot/to_plot.plg b/to_plot/to_plot.plg
new file mode 100755
index 0000000..14b9dc7
--- /dev/null
+++ b/to_plot/to_plot.plg
@@ -0,0 +1,24 @@
+
+
+
+Build Log
+
+--------------------Configuration: to_plot - Win32 Release--------------------
+
+Command Lines
+Creating temporary file "C:\DOCUME~1\NIEUWE~1\LOCALS~1\Temp\RSP3.tmp" with contents
+[
+/compile_only /nologo /warn:nofileopt /module:"Release/" /object:"Release/"
+"N:\simulations\dynamics.f90"
+]
+Output Window
+Compiling Fortran...
+N:\simulations\dynamics.f90
+
+
+
+Results
+dynamics.obj - 0 error(s), 0 warning(s)
+
+
+