      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='../../test_data/trapping/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
      CL=(PI4/DSQRT(X))*UIM**L
      MV=L+LF
      M=-LF
      DO 20 MF=1,MV
      M=M+1
      K=LFTL+M
      X=LFTL-M*(M+1)
      CP=DSQRT(X)
      X=LFTL-M*(M-1)
      CM=DSQRT(X)
      CZ=M
      W(K,1)=DCONJG(CP1*CP*YLM(K+2)+CM1*CM*YLM(K)+CZ1*CZ*YLM(K+1))*CL
   20 W(K,2)=DCONJG(CP2*CP*YLM(K+2)+CM2*CM*YLM(K)+CZ2*CZ*YLM(K+1))*CL
      DO 30 K=1,NLWM
      I=K+NLWM
      W(I,1)=UIM*W(K,2)
   30 W(I,2)=-UIM*W(K,1)
      RETURN
      END
CCC
