User talk:Monkeyjmp

fortran code for the force between two dielectric spheres

subroutine 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

for the 2d diffusion equation $u_{t}=u_{\text{xx}}+u_{\text{yy}}$

we have the following finite difference scheme

$\left(1-{\frac {1}{2}}\mu _{x}\delta _{x}^{2}\right)\left(1-{\frac {1}{2}}\mu _{y}\delta _{y}^{2}\right)U^{n+1}=\left(1+{\frac {1}{2}}\mu _{x}\delta _{x}^{2}\right)\left(1+{\frac {1}{2}}\mu _{y}\delta _{y}^{2}\right)U^{n}$

this is the peaceman rachford scheme it can be split into two schemes giving the adi (alternating direction implicit) form

$\left(1-{\frac {1}{2}}\mu _{x}\delta _{x}^{2}\right)U^{n+{\frac {1}{2}}}=\left(1+{\frac {1}{2}}\mu _{y}\delta _{y}^{2}\right)U^{n}$

$\left(1-{\frac {1}{2}}\mu _{y}\delta _{y}^{2}\right)U^{n}=\left(1+{\frac {1}{2}}\mu _{x}\delta _{x}^{2}\right)U^{n+{\frac {1}{2}}}$

images

[1] [2] [3] [4] [5] [6]

code for rp.f95

program 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