Upload New File

      PROGRAM FRFME

CCC   130331

CCC   FOCAL REGION FIELD (RICHARDS AND WOLF 1959);

CCC   INCIDENT GAUSSIAN PARAXIAL FIELD,

CCC   ABERRATION INDUCING SEPARATION PLANE ALLOWED WITH z = -SPD < 0

CCC   (NOVOTNY AND HECHT, NANO-OPTICS 2006)

CCC   LMODE=0 --> (0,0) MODE,

CCC   LMODE=1 -->  LPD  MODE,

CCC   LMODE=2 -->  RPD  MODE,

CCC   LMODE=3 -->  APD  MODE

      IMPLICIT REAL*8(A-H,O-Z)

CCC

CCC   GETS INPUT FOR BUILDING VK FROM TEDF BY EDFB OR FROM TMDF BY MDFB

CCC

      CHARACTER*4 NAMEF

      CHARACTER MORE

CCC   NLMMT=LM*(LM+2)*2;

CCC   NKV=NKS+1, NKS=NKSH*2, NRVC=NXV*NYV*NZV;

CCC   LM <= 20;

CCC   NKSH <= 500, NXSH, NYSH, NZSH <= 25

CCC   DIMENSION XV(NXV),YV(NYV),ZV(NZV),VKV(NKV),VKZM(NKV,NKV),

CCC  1WK(NLMMT),W(NKV,NKV),WSUM(NLMMT,NRVC)

      DIMENSION XV(51),YV(51),ZV(51),VKV(1001),VKZM(1001,1001),

     1WK(880),W(1001,1001),WSUM(880,132651)

      COMPLEX*16 WK,W,WSUM

      COMPLEX*16 CC0,SUMY,PHASF,PHASL,PHAS,SUMX

      IR=5

      IW=6

      IT=7

      ITT1=11

      ITT2=12

      CC0=(0.0D0,0.0D0)

      OPEN(IR,FILE='DFRFME',STATUS='OLD')

      READ(IR,*)JLMF,JLML

      IF(JLMF.EQ.1)GO TO 16

      CLOSE(IR)

      OPEN(IT,FILE='TFRFME',FORM='UNFORMATTED',STATUS='OLD')

      READ(IT)LMODE,LM,NKV,NXV,NYV,NZV

      READ(IT)VK,EXRI,AN,FF,TRA

      READ(IT)SPD,FRSH,EXRIL

      READ(IT)(XV(I),I=1,NXV)

      READ(IT)(YV(I),I=1,NYV)

      READ(IT)(ZV(I),I=1,NZV)

      OPEN(ITT2,FILE='TEMPTAPE2',FORM='UNFORMATTED',STATUS='OLD')

      READ(ITT2)(VKV(JX),JX=1,NKV)

      DO 10 JY=1,NKV

      DO 10 JX=1,NKV

      READ(ITT2)VKZM(JX,JY)

   10 CONTINUE

      READ(ITT2)APFAFA,PMF,SPD,RIR,FTCN,FSHMX,

     1VXYZMX,DELXYZ,VKNMX,DELK,DELKS,NLMMT,NRVC

      CLOSE(ITT2)

      DO 12 IXYZ=1,NRVC

      READ(IT)(WSUM(J,IXYZ),J=1,JLMF-1)

   12 CONTINUE

      CLOSE(IT)

      NKS=NKV-1

      GO TO 45

   16 READ(IR,*)LMODE,LM,NKSH,NRSH,NXSH,NYSH,NZSH,WLENFR

      READ(IR,*)AN,FF,TRA

      READ(IR,*)SPD,SPDFR,EXDCL

      READ(IR,*)MORE,IXI

      CLOSE(IR)

      IF(MORE.EQ.'m')MORE='M'

      IF(MORE.EQ.'e')MORE='E'

      IF((MORE.NE.'M').AND.(MORE.NE.'E'))GO TO 98

      NAMEF='T'//MORE//'DF'

      OPEN(IT,FILE=NAMEF,FORM='UNFORMATTED',STATUS='OLD')

      READ(IT)IDUML

      READ(IT)(IDUM,I=1,IDUML)

      READ(IT)EXDC,WP,XIP,IDFC,NXI

      IF(IDFC.LT.0)GO TO 18

      IF(IXI.GT.NXI)GO TO 96

      READ(IT)(XI,I=1,IXI)

      GO TO 20

   18 XI=XIP

   20 CLOSE(IT)

      WN=WP/3.0D08

      VK=XI*WN

      EXRI=DSQRT(EXDC)

      C0=0.0D0

      FRSH=C0

      EXRIL=C0

      FSHMX=C0

      APFAFA=EXRI/(AN*FF)

      IF(LMODE.NE.0)PMF=2.0D0*APFAFA

      IF(SPD.LE.C0)GO TO 22

      EXRIL=DSQRT(EXDCL)

      RIR=EXRI/EXRIL

      FTCN=2.0D0/(1.0D0+RIR)

      FRSH=-SPD*SPDFR

      STHMX=AN/EXRI

      STHLMX=STHMX*RIR

      UY=1.0D0

      FSHMX=SPD*

     1(RIR*(DSQRT(UY-STHMX*STHMX)/DSQRT(UY-STHLMX*STHLMX))-UY)

   22 NLMMT=LM*(LM+2)*2

      NKS=NKSH*2

      NKV=NKS+1      

      VKM=VK*EXRI

      VKNMX=VK*AN

      DELK=VKNMX/NKSH

      DELKS=DELK/VKM

      DELKS=DELKS*DELKS

      VXYZMX=(DACOS(C0)*4.0D0)/VKM

      VXYZMX=VXYZMX*WLENFR

      DELXYZ=VXYZMX/NRSH

      NXS=NXSH*2

      NXV=NXS+1

      NXSHPO=NXSH+1

      XV(NXSHPO)=C0

      DO 24 I=NXSHPO,NXS

      IPO=I+1

      XV(IPO)=XV(I)+DELXYZ

      XV(NXV-I)=-XV(IPO)

   24 CONTINUE

      NYS=NYSH*2

      NYV=NYS+1

      NYSHPO=NYSH+1

      YV(NYSHPO)=C0

      DO 25 I=NYSHPO,NYS

      IPO=I+1

      YV(IPO)=YV(I)+DELXYZ

      YV(NYV-I)=-YV(IPO)

   25 CONTINUE

      NZS=NZSH*2

      NZV=NZS+1

      NZSHPO=NZSH+1

      ZV(NZSHPO)=C0

      DO 27 I=NZSHPO,NZS

      IPO=I+1

      ZV(IPO)=ZV(I)+DELXYZ

      ZV(NZV-I)=-ZV(IPO)

   27 CONTINUE

      NRVC=NXV*NYV*NZV

      NKSHPO=NKSH+1

      VKV(NKSHPO)=C0

      DO 28 I=NKSHPO,NKS

      IPO=I+1

      VKV(IPO)=VKV(I)+DELK

      VKV(NKV-I)=-VKV(IPO)

   28 CONTINUE

      OPEN(IT,FILE='TFRFME',FORM='UNFORMATTED',STATUS='UNKNOWN')

      WRITE(IT)LMODE,LM,NKV,NXV,NYV,NZV

      WRITE(IT)VK,EXRI,AN,FF,TRA

      WRITE(IT)SPD,FRSH,EXRIL

      WRITE(IT)(XV(I),I=1,NXV)

      WRITE(IT)(YV(I),I=1,NYV)

      WRITE(IT)(ZV(I),I=1,NZV)

      OPEN(ITT1,FILE='TEMPTAPE1',FORM='UNFORMATTED',STATUS='UNKNOWN')

      OPEN(ITT2,FILE='TEMPTAPE2',FORM='UNFORMATTED',STATUS='UNKNOWN')

      WRITE(ITT2)(VKV(JX),JX=1,NKV)

      CALL FRFMER(NKV,VKM,VKV,VKNMX,APFAFA,TRA,

     1SPD,RIR,FTCN,LM,LMODE,PMF,ITT1,ITT2)

      CLOSE(ITT1)

      WRITE(ITT2)APFAFA,PMF,SPD,RIR,FTCN,FSHMX,

     1VXYZMX,DELXYZ,VKNMX,DELK,DELKS,NLMMT,NRVC

      CLOSE(ITT2)

      OPEN(ITT2,FILE='TEMPTAPE2',FORM='UNFORMATTED',STATUS='OLD')

      READ(ITT2)(VKV(JX),JX=1,NKV)

      DO 40 JY=1,NKV

      DO 40 JX=1,NKV

      READ(ITT2)VKZM(JX,JY)

   40 CONTINUE

      CLOSE(ITT2)

   45 DO 80 J=JLMF,JLML

      OPEN(ITT1,FILE='TEMPTAPE1',FORM='UNFORMATTED',STATUS='OLD')

      DO 50 JY=1,NKV

      DO 50 JX=1,NKV

      READ(ITT1)(WK(I),I=1,NLMMT)

      W(JX,JY)=WK(J)

   50 CONTINUE

      CLOSE(ITT1)

      IXYZ=0

      DO 75 IZ=1,NZV

      Z=ZV(IZ)+FRSH

      DO 70 IY=1,NYV

      Y=YV(IY)

      DO 65 IX=1,NXV

      X=XV(IX)

      IXYZ=IXYZ+1

      SUMY=CC0

      DO 60 JY=1,NKV

      VKY=VKV(JY)

      VKX=VKV(NKV)

      VKZF=VKZM(1,JY)

      PHASF=(0.0D0,1.0D0)*(-VKX*X+VKY*Y+VKZF*Z)

      PHASF=CDEXP(PHASF)

      VKZL=VKZM(NKV,JY)

      PHASL=(0.0D0,1.0D0)*(VKX*X+VKY*Y+VKZL*Z)

      PHASL=CDEXP(PHASL)

      SUMX=0.5D0*(W(1,JY)*PHASF+W(NKV,JY)*PHASL)

      DO 55 JX=2,NKS

      VKX=VKV(JX)

      VKZ=VKZM(JX,JY)

      PHAS=(0.0D0,1.0D0)*(VKX*X+VKY*Y+VKZ*Z)

      PHAS=CDEXP(PHAS)

      SUMX=SUMX+W(JX,JY)*PHAS

   55 CONTINUE

      IF((JY.EQ.1).OR.(JY.EQ.NKV))SUMX=0.5D0*SUMX

      SUMY=SUMY+SUMX

   60 CONTINUE

      WSUM(J,IXYZ)=SUMY*DELKS

   65 CONTINUE

   70 CONTINUE

   75 CONTINUE

   80 CONTINUE

      IF(JLMF.EQ.1)GO TO 88

      OPEN(IT,FILE='TFRFME',FORM='UNFORMATTED',STATUS='OLD')

      READ(IT)LMODE,LM,NKV,NXV,NYV,NZV

      READ(IT)VK,EXRI,AN,FF,TRA

      READ(IT)SPD,FRSH,EXRIL

      READ(IT)(XV(I),I=1,NXV)

      READ(IT)(YV(I),I=1,NYV)

      READ(IT)(ZV(I),I=1,NZV)

   88 DO 90 IXYZ=1,NRVC

      WRITE(IT)(WSUM(J,IXYZ),J=1,JLML)

   90 CONTINUE

      CLOSE(IT)

      OPEN(IW,FILE='OFRFME',STATUS='UNKNOWN')

      WRITE(IW,*)

     1'IF JLML < NLMMT, PRESERVE TEMPTAPE1, TEMPTAPE2, AND TFRFME,'

      WRITE(IW,*)'AND RESTART LM RUN WITH JLMF = JLML+1'

      IF(SPD.GT.0.0D0)WRITE(IW,*)' FSHMX =',FSHMX

      WRITE(IW,*)' FRSH =',FRSH

      CLOSE(IW)

      STOP

   96 CLOSE(IT)

   98 OPEN(IW,FILE='OFRFME',STATUS='UNKNOWN')

      WRITE(IW,*)' WRONG INPUT TAPE'

      CLOSE(IW)

      STOP

      END            

      SUBROUTINE FRFMER(NKV,VKM,VKV,VKNMX,APFAFA,TRA,

     1SPD,RIR,FTCN,LE,LMODE,PMF,ITT1,ITT2)

      IMPLICIT REAL*8(A-H,O-Z)

CCC   DIMENSION VKV(NKV)

      DIMENSION VKV(1001)

CCC   DIMENSION WK(NLEMT)

      DIMENSION WK(880)

      COMPLEX*16 WK

      COMPLEX*16 CC0

      CC0=(0.0D0,0.0D0)

      C0=0.0D0

      NLEMT=LE*(LE+2)*2

      SQ=VKM*VKM

      DO 90 JY=1,NKV

      VKY=VKV(JY)

      SQY=VKY*VKY

      DO 80 JX=1,NKV

      VKX=VKV(JX)

      SQX=VKX*VKX

      SQN=SQX+SQY

      VKN=DSQRT(SQN)

      IF(VKN.GT.VKNMX)GO TO 50

      VKZ=DSQRT(SQ-SQN)

      CALL WAMFF

     1(WK,VKX,VKY,VKZ,LE,APFAFA,TRA,SPD,RIR,FTCN,LMODE,PMF)

      WRITE(ITT1)(WK(J),J=1,NLEMT)

      WRITE(ITT2)VKZ

      GO TO 80

   50 WRITE(ITT1)(CC0,J=1,NLEMT)

      WRITE(ITT2)C0

   80 CONTINUE

   90 CONTINUE

      RETURN

      END

      SUBROUTINE WAMFF

     1(WK,X,Y,Z,LM,APFAFA,TRA,SPD,RIR,FTCN,LMODE,PMF)

      IMPLICIT REAL*8(A-H,O-Z)

CCC   DIMENSION WK(NLMMT)

CCC   NLMMT=NLMM*2, NLMM=LM*(LM+2)

      DIMENSION WK(880)

      COMPLEX*16 WK

CCC   DIMENSION W(NLMMT,2),YLM(NLMP)

CCC   NLMP=NLMM+2

      DIMENSION W(880,2),YLM(442),UP(3),UN(3)

      COMPLEX*16 W,YLM

      COMPLEX*16 CC0,UIM,CFAM,CF1,CF2

      LOGICAL ONX,ONY,ONZ

      CC0=(0.0D0,0.0D0)

      UIM=(0.0D0,1.0D0)

      ONE=1.0D0

      C0=0.0D0

      ONX=Y.EQ.C0

      ONY=X.EQ.C0

      ONZ=ONX.AND.ONY

      IF(ONZ)GO TO 10

      IF(ONX)GO TO 12

      IF(ONY)GO TO 13

      RHOS=X*X+Y*Y

      RHO=DSQRT(RHOS)

      CPH=X/RHO

      SPH=Y/RHO

      GO TO 15

   10 CPH=ONE

      SPH=C0

      GO TO 15

   12 RHOS=X*X

      RHO=DABS(X)

      CPH=DSIGN(ONE,X)

      SPH=C0

      GO TO 15

   13 RHOS=Y*Y

      RHO=DABS(Y)

      CPH=C0

      SPH=DSIGN(ONE,Y)

   15 IF(Z.NE.C0)GO TO 18

      IF(ONZ)GO TO 17

      R=RHO

      CTH=C0

      STH=ONE

      GO TO 22

   17 R=C0

      CTH=ONE

      STH=C0

      GO TO 22

   18 IF(ONZ)GO TO 20

      R=DSQRT(RHOS+Z*Z)

      CTH=Z/R

      STH=RHO/R

      GO TO 22

   20 R=DABS(Z)

      CTH=DSIGN(ONE,Z)

      STH=C0

   22 NLMM=LM*(LM+2)

      NLMMT=NLMM*2

      IF((LMODE.EQ.0).OR.(STH.NE.C0))GO TO 25

      DO 23 I=1,NLMMT

   23 WK(I)=CC0

      RETURN

   25 NLMP=NLMM+2

      YLM(NLMP)=(0.0D0,0.0D0)

      CALL SPHAR(CTH,STH,CPH,SPH,LM,YLM)

      UP(1)=CPH*CTH

      UP(2)=SPH*CTH

      UP(3)=-STH

      UN(1)=-SPH

      UN(2)=CPH

      UN(3)=C0

      CALL PWMALP(W,UP,UN,YLM,LM)

      APFA=STH*APFAFA

      IF(SPD.LE.C0)GO TO 57

      STHL=STH*RIR

      CTHL=DSQRT(1.0D0-STHL*STHL)

      ARG=SPD*(Z-(R/RIR)*CTHL)

      CFAM=(TRA*DEXP(-APFA*APFA)/DSQRT(CTHL))*CDEXP(UIM*ARG)

      GO TO (46,50,53),LMODE

      IF(STH.EQ.C0)GO TO 45

      GO TO 47

   45 FTC1=FTCN

      FTC2=FTCN

      GO TO 48

   46 CFAM=CFAM*(CPH+UIM*SPH)*(STH*PMF)

   47 FTC1=2.0D0*CTHL/(CTHL*RIR+CTH)

      FTC2=2.0D0*CTHL/(CTHL+CTH*RIR)

   48 CF1=(CPH*FTC1)*CFAM

      CF2=-(SPH*FTC2)*CFAM

      GO TO 62

   50 FTC1=2.0D0*CTHL/(CTHL*RIR+CTH)

      CFAM=CFAM*(STH*PMF*FTC1)

      DO 52 I=1,NLMMT

   52 WK(I)=W(I,1)*CFAM

      RETURN

   53 FTC2=2.0D0*CTHL/(CTHL+CTH*RIR)

      CFAM=CFAM*(STH*PMF*FTC2)

      DO 55 I=1,NLMMT

   55 WK(I)=W(I,2)*CFAM

      RETURN

   57 FAM=TRA*DEXP(-APFA*APFA)/DSQRT(CTH)

      GO TO (60,65,68),LMODE

      F1=CPH*FAM

      F2=-SPH*FAM

      DO 58 I=1,NLMMT

   58 WK(I)=W(I,1)*F1+W(I,2)*F2

      RETURN

   60 CFAM=(PMF*STH*FAM)*(CPH+UIM*SPH)

      CF1=CPH*CFAM

      CF2=-SPH*CFAM

   62 DO 63 I=1,NLMMT

   63 WK(I)=W(I,1)*CF1+W(I,2)*CF2

      RETURN

   65 FAM=PMF*STH*FAM

      DO 67 I=1,NLMMT

   67 WK(I)=W(I,1)*FAM

      RETURN

   68 FAM=PMF*STH*FAM

      DO 70 I=1,NLMMT

   70 WK(I)=W(I,2)*FAM

      RETURN

      END

      SUBROUTINE SPHAR(COSRTH,SINRTH,COSRPH,SINRPH,LL,YLM)

      IMPLICIT REAL*8(A-H,O-Z)

CCC   DIMENSION SINRMP(LL),COSRMP(LL),PLEGN((LLPO*LL)/2+LLPO),

CCC  1YLM(LLPO*LLPO)

      DIMENSION SINRMP(20),COSRMP(20),PLEGN(231),YLM(441)

      COMPLEX*16 YLM

      PI4=DACOS(0.0D0)*8.0D0

      PI4IRS=1.0D0/DSQRT(PI4)

      X=COSRTH

      Y=DABS(SINRTH)

      CLLMO=3.0D0

      CLL=1.5D0

      YTOL=Y

      PLEGN(1)=1.0D0

      PLEGN(2)=X*DSQRT(CLLMO)

      PLEGN(3)=YTOL*DSQRT(CLL)

      SINRMP(1)=SINRPH

      COSRMP(1)=COSRPH

      IF(LL.LT.2)GO TO 30

      K=3

      DO 20 L=2,LL

      LMO=L-1

      LTPO=L+L+1

      LTMO=LTPO-2

      LTS=LTPO*LTMO

      CN=LTS

      DO 10 MPO=1,LMO

      M=MPO-1

      MPOPK=MPO+K

      LS=(L+M)*(L-M)

      CD=LS

      CNM=LTPO*(LMO+M)*(L-MPO)

      CDM=(LTMO-2)*LS

   10 PLEGN(MPOPK)=PLEGN(MPOPK-L)*X*DSQRT(CN/CD)-

     1PLEGN(MPOPK-LTMO)*DSQRT(CNM/CDM)

      LPK=L+K

      CLTPO=LTPO

      PLEGN(LPK)=PLEGN(K)*X*DSQRT(CLTPO)

      K=LPK+1

      CLT=LTPO-1

      CLL=CLL*(CLTPO/CLT)

      YTOL=YTOL*Y

      PLEGN(K)=YTOL*DSQRT(CLL)

      SINRMP(L)=SINRPH*COSRMP(LMO)+COSRPH*SINRMP(LMO)

   20 COSRMP(L)=COSRPH*COSRMP(LMO)-SINRPH*SINRMP(LMO)

   30 L=0

   40 M=0

      K=L*(L+1)

      L0Y=K+1

      L0P=K/2+1

      YLM(L0Y)=PI4IRS*PLEGN(L0P)

      GO TO 45

   44 LMP=L0P+M

      SAVE=PI4IRS*PLEGN(LMP)

      LMY=L0Y+M

      YLM(LMY)=SAVE*DCMPLX(COSRMP(M),SINRMP(M))

      IF(MOD(M,2).NE.0)YLM(LMY)=-YLM(LMY)

      LMY=L0Y-M

      YLM(LMY)=SAVE*DCMPLX(COSRMP(M),-SINRMP(M))

 45   IF(M.GE.L)GO TO 47

      M=M+1

      GO TO 44

 47   IF(L.GE.LL) RETURN

      L=L+1

      GO TO 40

      END

      SUBROUTINE PWMALP(W,UP,UN,YLM,LW)

      IMPLICIT REAL*8(A-H,O-Z)

CCC   DIMENSION W(NLWMT,2),YLM(NLWM+2)

      DIMENSION W(880,2),UP(3),UN(3),YLM(442)

      COMPLEX*16 W,YLM

      COMPLEX*16 CP1,CM1,CP2,CM2,UIM,CL

      PI4=DACOS(0.0D0)*8.0D0

      NLWM=LW*(LW+2)

      UIM=(0.0D0,1.0D0)

      CM1=.5D0*DCMPLX(UP(1),UP(2))

      CP1=.5D0*DCMPLX(UP(1),-UP(2))

      CZ1=UP(3)

      CM2=.5D0*DCMPLX(UN(1),UN(2))

      CP2=.5D0*DCMPLX(UN(1),-UN(2))

      CZ2=UN(3)

      DO 20 L=1,LW

      LF=L+1

      LFTL=LF*L

      X=LFTL