! ! Strooiveld calculate the dipole field of an square array ! of dipoles. ! There are n * n dipoles, with spacing of d nm. The dipoles are assumed ! to have a thickness w. In other words: the volume is d² * w nm³ ! Dach dipole has a magnetic moment m0, with vector m. ! The size of m0 should be related to the volume via the density. ! ! The magnetic field will be calculated at a hight z above the film, ! and at random position x, y, where 0 < x,y < n * d. ! The weigh of te signal resulting from this field will be taken from ! the penetration calculation via TrimSP. ! ! ! $id: strooiveld.f90, v $ ! implicit none integer*4, parameter :: max_domain = 5000 ! sqaure root of maximum number of domains real*8 mr(3) ! position of the dipole. real*8 m0 ! magnitude of dipole real*8 domain(3,max_domain,max_domain) ! magnetization (vector) of domain real*8 f(3) ! position of field calculation real*8 b(3) ! resulting field real*8 gyro ! gyromagnetic cut of muon [MHz/tesla] real*8 factor ! factor to ensure field in tesla when the mommets ! are in Bohrmagneton and the spatial coordinates in ! nm. real*8 g_t(3,0:999) ! output depolarization curves ! (1,*) for forward-backward counters ! (2,*) for left-right counters ! (3,*) for top-bottom counters ! (*,0:999) from 0 to 9.99 microsecond real*8 initial(3) ! initial directioncosines of the muon real*8 d, w, domain_size ! size of the dipoles and the domains real*8 density ! number of magnetic atoms per cubic nm real*8 m_atom ! size the magnetic moment per atom integer*4 n ! number of dipolesalong one side real*8 angle ! angle between the dipole direction and x-axis integer*4 n_angle ! number of posible angles ! (e.g. if 1: ferromag, 2 up/down, 4 four-fold anisotropy) real*8 cut ! cut between the largest distance for which the ! dipole contribution should be calculated and the samllest ! distance. ! integer*4 i, j, k, lmuon, id, nd, ii, jj, imuon, nmuon, range, center_y, center_z real*8 r(3), r2, r3, r5, help, h(3), mx, my, mz, muon, square_field, current character*80 inputfile ! character*1 answer ! ! gyro = 135.5342D+00 ! MHz / tesla factor = 0.927412D-03 ! tesla * m³ / muBohr ! initial(1) = 0.0D+00 ! initial direction of the muons (x) initial(2) = 1.0D+00 ! initial direction of the muons (y) initial(3) = 0.0D+00 ! initial direction of the muons (z) cut = 1.0D+02 ! ! input.......................... (later to get from file or input) ! write(inputfile,123) 123 format('1000.rge') n = 5000 d = 10.0D+00 w = 20.0D+00 domain_size = 100.0D+00 n_angle = 360 density = 2.0D+00 / ( 0.28665D+00 ** 3 ) ! 2 atoms per bcc unit-cell of 2.866 Angstrom m_atom = 2.3D+00 ! assuming 2.3 Bohrmagneton per atom m0 = d * d * w * density * m_atom id = int( domain_size / d ) nd = n / id current = m0 * 9.27D-24 / ( d * d * 1.0D-18 ) write(6,*) 'm0 = ', m0, ' Bohrmagneton, magnetization / m² = ', current, ' Amp' write(6,*) 'If m=(1,0,0) then the field at center is ', square_field( current, float(n) * d, 0.0D+00 ) write(6,*) 'Press Enter to continue ' read(5,'(A1)') answer ! call plots() domain = 0.0D+00 DO i = 1, nd DO j = 1, nd ! angle = mod( int( float(n_angle) * rand() ), n_angle ) / float(n_angle) ! angle = 6.28D+00 * angle ! my = cos(angle) ! mz = sin(angle) ! mx = 1.0 ! my = 0.0 ! mz = 0.0 ! write(6,*) angle IF ( i .LT. nd / 4 ) THEN domain(2,i,j) = 1.0 domain(3,i,j) = 0.0 ELSE IF ( i .LT. nd / 2 ) THEN domain(2,i,j) = -1.0 domain(3,i,j) = 0.0 ELSE IF ( i .LT. 3 * nd / 4 ) THEN domain(2,i,j) = 0.0 domain(3,i,j) = -1.0 ELSE domain(2,i,j) = 0.0 domain(3,i,j) = 1.0 END IF END DO END DO domain = m0 * domain open(1, file=inputfile, status = 'old' ) DO i = 1, 7 read(1,'(A)') answer ! read first 7 lines with garbage END DO g_t = 0.0D+00 ! initialize depolarization nmuon = 0 1 read(1,*,end=2) f(1), muon ! read the x-position and number of muons stopped there lmuon = int( muon + 0.5) f(1) = f(1) / 1.0D+01 ! the .rge file from TrimSPL is in Angstrom IF ( f(1) .GT. w ) GOTO 2 ! stop is the muons go into the magnetic layer IF ( nmuon .EQ. 0 ) range = int( cut * ( 1.5D+00 * w - f(1) ) / d ) + 1 IF ( 2 * range .GT. n ) THEN Write(6,*) ' Sorry, range larger than n, have to stop ' STOP END IF write(6,*) 'Doing ', lmuon, ' muons at ', f(1), ' nanometer depth, with range of ', range lmuon = 1000 DO imuon = 1, lmuon ! loop over muons in one x-plane IF ( mod( imuon, 100 ) .EQ. 0 ) write(6,'(I5, 3F10.1,3D15.4)') imuon, f, b ! center_y = range + int(rand() * float(n - range - range) ) + 1 ! center_z = range + int(rand() * float(n - range - range) ) + 1 ! f(2) = (float(center_y)-0.5D+00) * d ! avoid to come too near to the edge ! f(3) = (float(center_z)-0.5D+00) * d ! of the sample f(2) = float(imuon) * d f(3) = 5.0D+03 ! center_y = int( f(2) / d + 0.5) ! center_z = int( f(3) / d + 0.5) b = 0.0D+00 DO i = 1, n DO j = 1, n ! DO i = center_y - range, center_y + range ! DO j = center_z - range, center_z + range mr(1) = 1.5D+00 * w ! the x-coordinate starts at the surface mr(2) = (float(i) - 0.5D+00) * d ! postion of center of magnetic moment mr(3) = (float(j) - 0.5D+00) * d r = mr - f ! vector from muon to magnetic moment r2 = sum( r * r ) ! square of length of r help = sqrt( r2) ! length of r r3 = help * r2 ! cube of r r5 = r2 * r3 ! r⁵ ii = 1 + i / id jj = 1 + j / id DO k = 1, 3 h(k) = domain(k,ii,jj) END DO help = sum( h * r ) ! scalar product of moment and r b = b + ( 3.0D+00 * help * r - r2 * h ) / r5 END DO END DO b = factor * b write(3,'(6D18.8)') f,b write(6,'(6D18.8)') f,b call precession( b, 0.0D+00, initial, g_t ) nmuon = nmuon + 1 END DO GOTO 1 2 close(1) g_t = g_t / float( nmuon ) ! normalize g_t open(2,file='muon.out',status='unknown') DO i = 0, 999 write(2,'(3F12.4)') (g_t(j,i),j=1,3) END DO close(2) ! np = np + 1 ! x(np) = sngl(d) * float(l) ! y1(np) = sngl(b(1)) ! y2(np) = sngl(b(2)) ! y3(np) = sngl(b(3)) ! END DO ! call plotxy(x, y1, np) ! call plotxy(x, y2, np) ! call plotxy(x, y3, np) ! ! write(6,3) f,b !3 format(' ', 6D12.2,' > '$) ! read(5,'(A1)') answer ! if ( answer .DQ. ' ') GOTO 1 STOP END !############################################################################################ SUBROUTINE precession( b, b_ext, emu, g_t ) real*8 b(3), emu(3), g_t(3,0:999), b_sq, b_abs, omega, t, ss, cc, eb(3), b_ext real*8 gyro, twpi integer*4 it data gyro/135.5D+00/ data twpi/6.283185307D+00/ ! 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.0D+00-cc)) * emu(1) + & & ( -eb(3)*ss+eb(1)*eb(2)*(1.0D+00-cc)) * emu(2) + & & ( eb(2)*ss+eb(1)*eb(3)*(1.0D+00-cc)) * emu(3) ! g_t(2,it) = g_t(2,it) + & & ( eb(3)*ss+eb(1)*eb(2)*(1.0D+00-cc)) * emu(1) + & & ( cc+eb(2)*eb(2)*(1.0D+00-cc)) * emu(2) + & & ( -eb(1)*ss+eb(2)*eb(3)*(1.0D+00-cc)) * emu(3) ! g_t(3,it) = g_t(3,it) + & & ( -eb(2)*ss+eb(1)*eb(3)*(1.0D+00-cc)) * emu(1) + & & ( eb(1)*ss+eb(2)*eb(3)*(1.0D+00-cc)) * emu(2) + & & ( cc+eb(3)*eb(3)*(1.0D+00-cc)) * emu(3) ! END DO RETURN END !############################################################################################# ! ! square_field returns the field of a square with sides l, carrying a current i at a heigth z ! above the center. ! current i in Ampere, l and z in nanometer ! real*8 FUNCTION square_field( i, l, z ) ! real*8 i, l, z, lm, zm ! ! convert to MKS ! lm = 1.0D-09 * l zm = 1.0D-09 * z ! square_field = 2.0D-07 * i * lm * lm / & & ((zm * zm + lm * lm / 4.0D+00) * sqrt(zm * zm + lm * lm / 2.0D+00)) RETURN END