electromag
editfortran code for the force between two dielectric spheres
editsubroutine spdivFsusee(R1,R2,k0,sdFee) implicit none !RETURNS !sdFee the spatial derivative of the !field suceptibility tensor for !the contribution to the electric field !at R1 due to dipole at R2 !sdFee is a 3x3x3 matrix !first index is the direction !of the derivative !next two indexes form a 3x3 !succeptibility tensor !for the spatial derivative !of the electric field !due to the electric field !of an electric dipole !FIXME !sdFee =(n, DxFxx DxFxy DxFxz ) ! ( DxFyx ... ... ! ( DxFxz ... ! ( DyFxx ... ! ( ... ... !================================= !ARGUMENTS !R1 is the current position !R2 is the position of the dipole !k0 is the scalar magnitude real*8, Dimension(3) :: R1,R2, R real*8 :: k0, Rmag, x,y,z complex*16, parameter :: ii=(0.,1.) complex*16 :: A,B,C,D complex*16, dimension(3,3,3), intent(out) :: sdFee R=R1-R2 x=R(1) y=R(2) z=R(3) Rmag= Sqrt(R(1)*R(1)+R(2)*R(2)+R(3)*R(3)) A=k0**2*((ii/(k0*2*Rmag**2))-(4./(k0**3*Rmag**3))-& (9.*ii/(k0**4*Rmag**4))+(9./(k0**5*Rmag**5))) B=k0**4*((-ii/(k0**4*Rmag**4))+(6./(k0**5*Rmag**5))+& (15.*ii/(k0**6*Rmag**6))-(15./(k0**7*Rmag**7))) C=k0**2*((ii/(k0**2*Rmag**2))-(2./(k0**3*Rmag**3))-& (3.*ii/(k0**4*Rmag**4))+(3./(k0**5*Rmag**5))) D=k0**2*((-1./(k0**3*Rmag**4))-(3.*ii/(k0**4*Rmag**4))+& (3./(k0**5*Rmag**5))) !write (*,*) 'A B C D', A, B, C, D sdFee=0. sdFee(1,1,1)=A*x+B*x**3 sdFee(1,1,2)=D*y+B*x**2*y sdFee(1,1,3)=D*z+B*x**2*z sdFee(1,2,1)=sdFee(1,1,2) sdFee(1,2,2)=C*x+B*x*y**2 sdFee(1,2,3)=B*x*y*z sdFee(1,3,1)=sdFee(1,1,3) sdFee(1,3,2)=sdFee(1,2,3) sdFee(1,3,3)=C*x+B*x*z**2 !================================== sdFee(2,1,1)=C*y+B*x**2*y sdFee(2,1,2)=D*x+B*x*y**2 sdFee(2,1,3)=B*y*x*z sdFee(2,2,1)=sdFee(2,1,2) sdFee(2,2,2)=A*y+B*y**3 sdFee(2,2,3)=D*z+B*y**2*z sdFee(2,3,1)=sdFee(2,1,3) sdFee(2,3,2)=sdFee(2,2,3) sdFee(2,3,3)=C*y+B*y*z**2 !================================ sdFee(3,1,1)=C*z+B*x**2*z sdFee(3,1,2)=B*z*x*y sdFee(3,1,3)=D*x+B*x*z**2 sdFee(3,2,1)=sdFee(3,1,2) sdFee(3,2,2)=C*z+B*y**2*z sdFee(3,2,3)=D*y+B*y*z**2 sdFee(3,3,1)=sdFee(3,1,3) sdFee(3,3,2)=sdFee(3,2,3) sdFee(3,3,3)=A*z+B*z**3 sdFee=k0**3 * exp(ii*k0*Rmag)*sdFee ! write (*,*) 'sdFee',sdFee(1,1,1), sdFee(1,2,2), sdFee(1,3,3) end subroutine subroutine E_iwave(R,E0,k0,E_i) implicit none !RETURNS !E_i is the electromagnetic field !due to the incident plane wave !================================== !ARGUMENTS !k0 is the wave vector !should be perpendicular to E0 !for the incident plane wave !E(r) = E_0 exp(i(k . r)) real*8 , dimension(3) :: R,E0, k0,tmp complex*16, dimension(3), intent(out) :: E_i complex*16, parameter :: ii=(0.,1.) tmp=dot_product(k0,R) E_i = E0 * exp(ii*tmp) end subroutine E_iwave real function Rmagnitude(R) implicit none real*8 , dimension(3) :: R Rmagnitude = Sqrt(R(1)*R(1)+R(2)*R(2)+R(3)*R(3)) end function Rmagnitude subroutine FSus(R1,R2,k0,Fee) implicit none !RETURNS !Fee the field suceptibility tensor !at R1 due to dipole at R2 !================================= !ARGUMENTS !R1 is the current position !R2 is the position of the dipole !k0 is the scalar magnitude real*8, Dimension(3) :: R1,R2, R real*8 :: k0, Rmag integer :: i,j, kdelta complex*16, parameter :: ii=(0.,1.) complex*16, dimension(3,3), intent(out) :: Fee R=R1-R2 Rmag=Sqrt(R(1)*R(1)+R(2)*R(2)+R(3)*R(3)) Fee=0. do i=1,3 do j=1,3 Fee(i,j)=k0**3 * exp(ii*k0*Rmag)*& ((1./(k0*Rmag)+ii/(k0**2*Rmag**2)-& 1./(k0**3.*Rmag**3))*kdelta(i,j)-& (1./(k0*Rmag)+3.*ii/(k0**2*Rmag**2)& -3./(k0**3*Rmag**3))*(R(i)*R(j)/Rmag**2)) enddo enddo end subroutine FSus integer function kdelta(a,b) implicit none !the kronecker delta function integer a,b if (a==b) kdelta=1 if (a/=b) kdelta=0 end function kdelta !!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!! !BEGINING OF PROGRAM! !!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!! program force2 implicit none integer, parameter :: n=2 real*8, dimension(3,n) :: R_n complex*16, dimension(3,n) :: E_n real*8, dimension(3) :: kvec0, E0, force complex*16, dimension(3*n) :: E_inc,E_res complex*16, dimension(3*n,3*n) :: Gee, invGee complex*16, dimension(3,3) :: Idmat complex*16, dimension(3,3,n,n) :: Fee_n complex*16, dimension(3,3,3) :: sdFee !E_inc is the vector which holds all the E_incident's st !G*E_res = E_inc !E_res is the vector to hold the E_resultant's st !E_res = G^-1*E_inc integer i,j,kdelta,s,t,u,v,spaceloop, maxspaceloop complex*16, dimension(n) :: alpha_n complex*16 :: ii, alpha, alpha0 real*8 :: pi, k0, spacegran, radius, permitivity1, permitivity2 !needed for the code to walk the plane !and print out the values of the field ! real*8 :: x,y,z,startx,endx,starty,endy !=========================== !Vars for inversion ONLY integer :: info integer, dimension(3*n) :: ipivot complex*16, dimension(3*n) :: work !=========================== pi=acos(-1.) ii=(0.,1.) !for a transverse plane !wave propagating in the !x direction and !polarised in the z axis !k0 is wave number= 2pi/wavelength !kvec0 is the wave vector E0=(/0.d0,1.d0,0.d0/) kvec0=(/0.d0,0.d0,1.d0/) !for a wavelength of 632.8nm k0=2.*pi*sqrt(1.69)/(632.8d-9) !k0=2.*pi/1. !(632.8*(10.**(-9))) kvec0=k0*kvec0 maxspaceloop=5000 !write (*,'(I6)') maxspaceloop !sphere B is sphere number 1 !and is at R_n(:,1) !the position is sphere B is varied !in the graph to be reproduced spacegran=0.01d-9 R_n(:,1)=(/0.d0,0.d0,0.d0/) R_n(:,2)=(/0.,0.,0./) !alpha = alpha0 / (1-2/3 i k0^3 alpha0) !alpha0 = a^3 (e-1)/(e+2) !radius = 10nm radius = 10.d-9 permitivity1=2.25d0 /1.69 permitivity2=2.25d0 /1.69 alpha0 = (radius**3) * (permitivity1 -1.)/(permitivity1+2.) alpha=alpha0/(1.- (2./3.)*ii*(k0**3)*alpha0) alpha_n(1)=alpha ! alpha_n(1)=0.04 ! alpha_n(2)=0.04 !write (*,*) alpha, k0 alpha0 = (radius**3) * (permitivity2 -1.)/(permitivity2+2.) alpha=alpha0/(1.- (2./3.)*ii*(k0**3)*alpha0) alpha_n(2)=alpha ! write (*,*) alpha_n(1), alpha_n(2) do spaceloop=1,maxspaceloop R_n(1,1)=100.d-9+spaceloop*(2000.d-9-100.d-9)/REAL(maxspaceloop) !construct the Gee after constructing Fee_n !Fee_n(.,.,i,j) at R_n(i) due to dipole p_j=alpha_j*En_j at Rn_j !for each i,j = 1..n Fee_n=0. do i=1,n do j=1,n if (kdelta(i,j)==0) then Call FSus(R_n(:,i),R_n(:,j),k0,Fee_n(:,:,i,j)) endif enddo enddo !We construct Gee, for n=2 ! Gee = (I Fee_12 ) ! (Fee_21 I ) Idmat=0. do i=1,3 Idmat(i,i)=(1.,0.) enddo do i=1,n do j=1,n s=3*i-2 t=3*i u=3*j-2 v=3*j if (kdelta(i,j)==0) then Gee(s:t,u:v)=-1.*Fee_n(:,:,i,j)*alpha_n(j) else Gee(s:t,u:v)=Idmat endif enddo enddo !Construct E_inc the vector !which holds all the E_incident's do i=1,n s=3*i-2 t=3*i call E_iwave(R_n(:,i),E0,kvec0,E_inc(s:t)) enddo !NEED matrix inversion, !Gee -> invGee !trying ZGETRI and ZGETRF !from lapack 3.1.1 invGee=Gee call ZGETRF(3*n, 3*n, invGee, 3*n, IPIVOT, INFO) call ZGETRI(3*n, invGee, 3*n, IPIVOT, work, 3*n, INFO) !Calculate the resultant E field !E_res = G^-1*E_inc !these are all complex valued E_res=matmul(invGee,E_inc) !unpack the E_res into E_n do i=1,n s=3*i-2 t=3*i E_n(:,i)=E_res(s:t) enddo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !Walk the x-y plane and return tripples ! x,y,E_tot ! where E is determined by ! = pos = (x,y) ! E_tot = Sum_n F_pos(i,j)|i>j} p_i !u=200 !v=200 !startx=-18. !endx=18. !starty=-18. !endy=18. !!write the size of the grid !write (*, '(2X,2(I6,1X))') u,v !z=0.5 !do i=1,u ! x=startx+Real(i-1)*(endx-startx)/Real(u-1) ! do j=1,v ! y=starty+Real(j-1)*(endy-starty)/Real(v-1) ! Etmp=0. ! Pos = x*(/1.,0.,0./)+y*(/0.,1.,0./)+z*(/0.,0.,1./) ! do s=1,n ! Call FSus(Pos,R_n(:,s),k0,Ftmp) ! Etmp = Etmp + matmul(Ftmp,E_res(3*s-2:3*s)) ! !write (*,*) "x y etmp" ! ! enddo ! call E_iwave(Pos,E0,kvec0,Etmp2) ! Etmp=Etmp + Etmp2 ! Write (*,'(2X,(3E12.4,1X))') x,y,(Real(Etmp(1))) ! enddo !enddo !!!!!!!!!!!!!!!!!!!!!!!!! !These formatting strings !are special cased for n=2 only !write (*,*) "This is Gee" !write (*, '(6(E10.2,1X,E10.2,2X))') Gee !write (*,*) !write (*,*) "InvGee" !write (*, '(6(E10.2,1X,E10.2,2X))') InvGee !write (*,*) !write (*,*) "E0 E_inc E_res" !write (*,'(6(E9.1,1X)/)') E_inc !write (*,'(6(E9.1,1X)/)') E_res !!!!!!!!!!!!!!!!!!!!!!!!!! !Calculate the forces for the !using !<F^i> = 1/2 Re[alpha_1 E_1j partial^i(E_0^j)*] !<F^x> = 1/2 Re[alpha_1 E_1j alpha_2* E_2j* partial_x T12*] !for the x-force we require !the x-derivative of sdFee force=0. !t loops over the components of the force ! do t=1,3 t=3 s=2 ! Call FSus(R_n(:,1),R_n(:,2),k0,Fee_n(:,:,1,2)) call spdivFsusee(R_n(:,1),R_n(:,2),k0,sdFee) ! force(1)=0.5*real(alpha_n(2)*CONJG(alpha_n(1))*& ! E_n(s,2)*CONJG(E_n(s,1))*CONJG(sdFee(1,s,s))) force(1)=0.5*real(alpha_n(2)*CONJG(alpha_n(1))*& E_n(s,2)*CONJG(E_n(s,1))*CONJG(sdFee(1,s,s)) + & alpha_n(2)*CONJG(alpha_n(1))*E_n(s,2)*CONJG(ii*kvec0(1)*E_n(s,1))) ! enddo write (20,'(2E15.4)') R_n(1,1), force(1) enddo end program force2
finite difference stuff
editfor the 2d diffusion equation
we have the following finite difference scheme
this is the peaceman rachford scheme it can be split into two schemes giving the adi (alternating direction implicit) form
images
editcode for rp.f95
editprogram rp implicit none integer stepsj,stepsn integer n integer i,j real*8, dimension(:,:), allocatable :: w_now, w_next, w_half real*8, dimension(:,:) , allocatable :: w_row1,w_row2 !as deltax=deltay, we only need only need one mu real*8 :: mu !w_row1,w_row2 contain for each j=0..J !(:,1)=a_j !(:,2)=b_j !(:,3)=c_j !(:,4)=d_j !(:,5)=e_j !(:,6)=f_j ! stepsj=10 ! stepsn=10 open(unit=1,file='ic.pgm') open(unit=2,file='control.cfg') read(unit=1,'(/,/,I3,1X,I3)') stepsj,stepsj read (unit=2,*) stepsn,mu stepsj=stepsj+1 ! mu=1./(real(stepsj)**2) allocate(w_now(0:stepsj,0:stepsj)) allocate(w_next(0:stepsj,0:stepsj)) allocate(w_half(0:stepsj,0:stepsj)) allocate(w_row1(0:stepsj,6)) allocate(w_row2(0:stepsj,6)) !first we populate the w_now with the IC !this represents the u(x,y,0) !set to zero clears mem and !also sets the dirichlet BC at the edges w_now=0. read (1,*) w_now(1:stepsj-1,1:stepsj-1) !time loop do n=0,stepsn-1 !after all stepsn time steps !we output the w_now write (30+n,'(I8)') stepsj write (30+n,'(I8)') stepsj do i=0,stepsj do j=0,stepsj write (30+n,'(1X,3(E12.5,X))') real(i)/real(stepsj),real(j)/real(stepsj),w_now(i,j) enddo enddo !set to zero clears mem and sets BC's w_row1=0. w_half=0. do j=1,stepsj-1 do i=1,stepsj-1 !populate w_row1 with a_i,b_i,c_i,d_i w_row1(i,1)=.5*mu w_row1(i,2)=1.+mu w_row1(i,3)=.5*mu w_row1(i,4)=.5*mu*w_now(i,j-1)+& (1.-.5*mu)*w_now(i,j)+& .5*mu*w_now(i,j+1) !calculate e_i and f_i !e_i=c_i/(b_i-a_i*e_(i-1)) !f_i=(d_i+a_i*f_(i-1))/(b_i-a_i*e_(i-1)) w_row1(i,5)=w_row1(i,3)& /(w_row1(i,2)-w_row1(i,1)*w_row1(i-1,5)) w_row1(i,6)=(w_row1(i,4)+w_row1(i,1)*w_row1(i-1,6))& /(w_row1(i,2)-w_row1(i,1)*w_row1(i-1,5)) enddo !using the a_i..f_i we !calculate the jth row of the half time step do i=stepsj-1,1,-1 !U_i=f_i+e_i*U_(i+1) w_half(i,j)=w_row1(i,6)+w_row1(i,5)*w_half(i+1,j) enddo enddo !2nd phase, calculate w_next from w_half !set to zero clears mem and sets BC's w_row2=0. w_next=0. do i=1,stepsj-1 do j=1,stepsj-1 !populate w_row2 with a_j,b_j,c_j,d_j w_row2(j,1)=.5*mu w_row2(j,2)=1.+mu w_row2(j,3)=.5*mu w_row2(j,4)=.5*mu*w_half(i-1,j)+& (1.-.5*mu)*w_half(i,j)+& .5*mu*w_half(i+1,j) !calculate e_j and f_j !e_j=c_j/(b_j-a_j*e_(j-1)) !f_j=(d_j+a_j*f_(j-1))/(b_j-a_j*e_(j-1)) w_row2(j,5)=w_row2(j,3)& /(w_row2(j,2)-w_row2(j,1)*w_row2(j-1,5)) w_row2(j,6)=(w_row2(j,4)+w_row2(j,1)*w_row2(j-1,6))& /(w_row2(j,2)-w_row2(j,1)*w_row2(j-1,5)) enddo !using the a_j..f_j we !calculate the ith column of the next (or full) time step do j=stepsj-1,1,-1 !U_j=f_j+e_j*U_(j+1) w_next(i,j)=w_row2(j,6)+w_row2(j,5)*w_next(i+1,j) enddo enddo w_now=w_next enddo end program rp