From af60175bd6b2aaceaf7b1ca5cc1e73fed1d9a504 Mon Sep 17 00:00:00 2001 From: Giacomo Mulas <giacomo.mulas@inaf.it> Date: Mon, 4 Sep 2023 12:21:14 +0000 Subject: [PATCH] Upload New File --- trapping/frfme.f | 496 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 496 insertions(+) create mode 100644 trapping/frfme.f diff --git a/trapping/frfme.f b/trapping/frfme.f new file mode 100644 index 00000000..3c368a4a --- /dev/null +++ b/trapping/frfme.f @@ -0,0 +1,496 @@ + 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 + 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 \ No newline at end of file -- GitLab