!  dynamics.f90 
!
!
!****************************************************************************
!
!  PROGRAM: dynamics
!
!  PURPOSE: Simulation of the asymmetry of an artificial spinglass.
!			DYNAMICS assumes the spinglass to have an FCC lattice. 
!			The dimensions of the lattice are w * w * d, along the 
!			x-, y- and z-axis respectively. The magnetic moments
!			are randomly distributed over the latticepoints, the muons
!			are placed on the centers of the FCC-cube.
!			The directions of the magnetic moments is choosen randomly
!			over the whole sphere.
!			The program calculates the magnetic field at the site of 
!			muon by adding all dilopar contributions from about 300 
!			magnetic moments which nearest by the muonsite. Periodic
!			boundary conditions are applied in the x- and y-direction. 
!			The z-direction is assumed to perpendicular to a thin
!			film surface. 	
!			The dynamics of the magnetic spins is included in one of 
!			the following ways:	
!			For fluctuationrates larger then 100 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 10 microsecods or more.
!			Output of the muon position is done about every time_resolution microsecond.
!			For fluctuationrates smaller then 100 MHz
!			the muons rotate 1000 times for time_resolution microsecond, after each rotation
!			a fraction (= fluctuationrate / 100) 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
!
!	USE:	The parameters used for the simulation are supposed to be on
!			file with the generic name <calculation>.inp.
!			The program can be started in two ways:
!
!				typing DYNAMICS 
!					the user will be prompted for the name of the calculation
!
!				typing DYNAMICS <calculation>
!					the name of the calculation will be read from the commandline.
!
!			Output will be written on <calculation>.out and on separate files
!			(for each set of parameters) named <calculation>_###.g_t, where
!			### can a unique number according to the following rules:
!			If a file \simulations\counter.his can be opened, the program will
!			the number in this file and uses that as a start for numbering
!			the *.g_t files. The program will update \simulations\counter.his.
!			If that file is not present, the program will start at number 1.
!
!	INPUT:	For each simulation the following set of parameters has to be
!			given on one line in the file <calculation>.inp
!
!			lattice parameter							[nm]
!			magnetic moment								[Bohr-magneton]
!			external field, three component				[tesla]
!			thickness d									[nm]
!			width w										[nm]
!			concentration								[at.%]
!			number of muons								#
!			initial muon spin direction in 
!				spherical coordinates, theta, phi		[degree]
!				note that the z-axis is perpendicular to the film
!			muon stopping range, from d1 to d2			[nm]
!				d1 and d2 are note restricted by 0 and d, stopping outside
!				the actual sample is possible
!			fluctuationrate								[ 1/ microsecond ]
!			fluctuation amplitude, in 
!				spherical coordinates, d-theta, d-phi	[degree]
!
!			Lines with parameters can be interlaced with comments, 
!			Commentlines should have a  !  at position 1. 
!			
!
!
!****************************************************************************

	program dynamics

	Use DFPORT
	Use DFLIB

	implicit none

	! Variables

	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 
 &						 n_time_steps	=   500     ! number of time steps calculated
!													  Should be future variable
!													  Should be future variable


!	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    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

!	Declarations, maximumnumber of spins: max_spins, maxd is the maximum number of
!   unitcell-distance for which the spin in included in the calculation

	real*8,		  time_resolution=0.01		! time resolution 
!											(approximate time between points
!											 on the *.g_t file)


	character*10  dddd, tttt, zone
	character*4   file_index
	integer*4     dt(8), ifile, l_calc, n_steps, i_step
	character*80  calculation
	character*127 line
	logical		  unique
	integer*4     i,j,k,l,nw,nd,nsp,n_spin, n_site
	integer*4     iseed, nd1, nd2
	record        /spin/  s(max_spins)
	record        /muon/  m(max_muons)
	real*8        dummy, a, d, concentration, w, depth1, depth2
	real*8        factor, moment, b_ext(3)
	real*8        fluctuationrate, tau, dphi, dtheta, fraction
	real*8        g_t(3), omega, b_abs, b_sq, ca_sq, his, radiussq
	real*8        t_ini, p_ini, emu(3), Pi
	real*8		  step, exp_time, write_time, f_c

	! Body of dynamics

	! Read the parameters from the input file 
	! with name : <calculation>.inp
	! The output will go to <calculation>.out
	! and <calculation>_###.g_t

	Write(6,*) ' '
	Write(6,*) ' ---------------------------------------------------------------------' 
	Write(6,*) ' | Program field-calculation of muons due to random dynamic spins    |'
	Write(6,*) ' |                                                                   |'
	Write(6,*) ' |                 Version of November 16 2005                       |'
	Write(6,*) ' |                                                                   |'
	Write(6,*) ' | Input is read from an input file that should be named             |'
	Write(6,*) ' |     <calculation>.inp     and contains for each simulation on     |'
	Write(6,*) ' |                           one line:                               |'
	Write(6,*) ' |                                                                   |'
	Write(6,*) ' |              lattice-constant [nm], magnetic moment [mu_B],       |'
	Write(6,*) ' |              ext. field(3) ,thickness, width, c, number_of_muons, |'
	Write(6,*) ' |              initial-muon-direction(theta, phi)[degrees],         |'
	Write(6,*) ' |              (muon-positions from) depth1, (to) depth2 [nm],      |'
	Write(6,*) ' |              fluctions rate [inverse microsec],                   |'
	Write(6,*) ' |              fluctuation amplitude parallel to film [0..360degr.],|'
	Write(6,*) ' |              fluctuation amplitude perpendicular to film          |'
	Write(6,*) ' |                                                     [0..180degr.] |'
	Write(6,*) ' |                                                                   |'
	Write(6,*) ' | Lines with a ! in the first position are treated as comments.     |'
	Write(6,*) ' |                                                                   |'
	Write(6,*) ' | <calculation> 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