From c9244503f69a5abcee49e1cbb18f6724b4f9161b Mon Sep 17 00:00:00 2001 From: Giacomo Mulas <giacomo.mulas@inaf.it> Date: Mon, 4 Sep 2023 12:21:32 +0000 Subject: [PATCH] Upload New File --- trapping/lffft.f | 494 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 494 insertions(+) create mode 100644 trapping/lffft.f diff --git a/trapping/lffft.f b/trapping/lffft.f new file mode 100644 index 00000000..a805621c --- /dev/null +++ b/trapping/lffft.f @@ -0,0 +1,494 @@ + PROGRAM LFFFT +CCC 140630 +CCC LE=20, LM=20, LE<=LM +CCC FRPIF IS INTENSITY FACTOR ANALOGOUS TO I0 FOR FOCAL REGION POINTS; +CCC JFT <= 0 --> OUTPUTS ON TAPE FORCE*(c/FRPIF) IN m**2; +CCC JFT >= 0 --> OUTPUTS ON TAPE TORQUE*(c/FRPIF)*(VK*EXRI)**3; +CCC (TORQUE OUTPUTS)/(EXDC*VK*VK) ARE IN m**2 +CCC JSS = 1 USES TWS FOR INPUT +CCC JTW = 1 FOR OUTPUTS ON PAPER + IMPLICIT REAL*8(A-H,O-Z) +CCC DIMENSION ZPV(LE,3,2,2),AC(NLEMT),WS(NLEMT) + DIMENSION ZPV(20,3,2,2),AC(880),WS(880), + 1FFFE(3),FFFS(3),FFTE(3),FFTS(3) + COMPLEX*16 AC,WS +CCC DIMENSION XV(NXV),YV(NYV),ZV(NZV),WSL(NLMMT) + DIMENSION XV(51),YV(51),ZV(51),WSL(160) + COMPLEX*16 WSL +CCC DIMENSION AM0M(NLEMT,NLEMT) + DIMENSION AM0M(880,880) + COMPLEX*16 AM0M +CCC DIMENSION AMD(NVAM,4),INDAM(LE,LE) +CCC NVAM=LE*LE+(LE*LEPO*LETPO)/3 + DIMENSION AMD(6140,4),INDAM(20,20) + COMPLEX*16 AMD +CCC DIMENSION TMSM(LE),TMSE(LE) + DIMENSION TMSM(20),TMSE(20) + COMPLEX*16 TMSM,TMSE +CCC DIMENSION TMS(LE,3) + DIMENSION TMS(20,3) + COMPLEX*16 TMS + COMMON/CIL/LE,NLEM,NLEMT,MXMPO,MXIM + COMMON/CCR/CC0,UIM,SQ2ITI,COF,CIMU,SQ2I + COMPLEX*16 CC0,UIM,SQ2ITI + IR=5 + IT=7 + OPEN(IR,FILE='DLFFFT',STATUS='OLD') + READ(IR,*)JFT,JSS,JTW + CLOSE(IR) + OPEN(IT,FILE='TTMS',FORM='UNFORMATTED',STATUS='OLD') + READ(IT)IS,LE + READ(IT)VKS,EXRIS + NLEM=LE*(LE+2) + NLEMT=NLEM+NLEM + IF(IS.GT.2222)GO TO 120 + IF(IS.GT.1111)GO TO 125 + IF(IS.GE.0)GO TO 135 + IF(IS.LT.0)GO TO 145 + 120 CONTINUE + READ(IT)(TMS(I,1),TMS(I,2),TMS(I,3),I=1,LM) + GO TO 150 + 125 CONTINUE + READ(IT)(TMSM(I),TMSE(I),I=1,LE) + GO TO 150 + 135 CONTINUE + READ(IT)((AM0M(I,J),J=1,NLEMT),I=1,NLEMT) + GO TO 150 + 145 CONTINUE + NVAM=LE*LE+(LE*(LE+1)*(LE*2+1))/3 + READ(IT)((AMD(I,J),J=1,4),I=1,NVAM) + READ(IT)((INDAM(I,J),J=1,LE),I=1,LE) + READ(IT)MXMPO + MXIM=MXMPO*2-1 + 150 CLOSE(IT) + IF(JSS.NE.1) + 1OPEN(IT,FILE='TFRFME',FORM='UNFORMATTED',STATUS='OLD') + IF(JSS.EQ.1) + 1OPEN(IT,FILE='TWS',FORM='UNFORMATTED',STATUS='OLD') + READ(IT)LMODE,LM,NKV,NXV,NYV,NZV + IF(LM.LT.LE)GO TO 500 + READ(IT)VK,EXRI,AN,FF,TRA + IF((VK.NE.VKS).OR.(EXRI.NE.EXRIS))GO TO 500 + 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) + CC0=(0.0D0,0.0D0) + UIM=(0.0D0,1.0D0) + IF(JFT.GT.0)GO TO 155 + CALL THDPS(LE,ZPV) + EXDC=EXRI*EXRI + SQK=VK*VK*EXDC + COF=1.0D0/SQK + CIMU=COF/DSQRT(2.0D0) + IF(JSS.EQ.1)GO TO 155 + ITF=8 + OPEN(ITF,FILE='TLFFF',FORM='UNFORMATTED',STATUS='UNKNOWN') + WRITE(ITF)LMODE,LE,NKV,NXV,NYV,NZV + WRITE(ITF)VK,EXRI,AN,FF,TRA + WRITE(ITF)SPD,FRSH,EXRIL + WRITE(ITF)(XV(I),I=1,NXV) + WRITE(ITF)(YV(I),I=1,NYV) + WRITE(ITF)(ZV(I),I=1,NZV) + IF(JFT.LT.0)GO TO 160 + 155 SQ2I=1.0D0/DSQRT(2.0D0) + SQ2ITI=SQ2I*UIM + IF(JSS.EQ.1)GO TO 160 + ITT=9 + OPEN(ITT,FILE='TLFFT',FORM='UNFORMATTED',STATUS='UNKNOWN') + WRITE(ITT)LMODE,LE,NKV,NXV,NYV,NZV + WRITE(ITT)VK,EXRI,AN,FF,TRA + WRITE(ITT)SPD,FRSH,EXRIL + WRITE(ITT)(XV(I),I=1,NXV) + WRITE(ITT)(YV(I),I=1,NYV) + WRITE(ITT)(ZV(I),I=1,NZV) + 160 NLMM=LM*(LM+2) + NLMMT=NLMM+NLMM + DO 475 IZ=1,NZV + DO 475 IY=1,NYV + DO 475 IX=1,NXV + IF(LM.GT.LE)GO TO 170 + READ(IT)(WS(I),I=1,NLMMT) + GO TO 180 + 170 READ(IT)(WSL(I),I=1,NLMMT) + DO 175 I=1,NLEM + IE=I+NLEM + IEL=I+NLMM + WS(I)=WSL(I) + WS(IE)=WSL(IEL) + 175 CONTINUE + 180 IF(IS.EQ.2222)GO TO 200 + IF(IS.EQ.1111)GO TO 205 + IF(IS.GT.0)GO TO 305 + IF(IS.LT.0)GO TO 405 + 200 CONTINUE + CALL SAMPOA(AC,TMS,WS) + GO TO 445 + 205 CONTINUE + CALL SAMP(AC,TMSM,TMSE,WS) + GO TO 445 + 305 CONTINUE + CALL CAMP(AC,AM0M,WS) + GO TO 445 + 405 CONTINUE + CALL CZAMP(AC,AMD,INDAM,WS) + 445 CONTINUE + IF(JFT.GT.0)GO TO 460 + CALL FFRF(ZPV,AC,WS,FFFE,FFFS) + IF(JSS.NE.1)GO TO 450 + WRITE(66,*)FFFE(1),FFFS(1),FFFE(1)-FFFS(1) + WRITE(66,*)FFFE(2),FFFS(2),FFFE(2)-FFFS(2) + WRITE(66,*)FFFE(3),FFFS(3),FFFE(3)-FFFS(3) + GO TO 455 + 450 WRITE(ITF)(FFFE(I)-FFFS(I),I=1,3) + IF(JTW.EQ.1)WRITE(66,6666) + 1IX,IY,IZ,FFFE(1)-FFFS(1),FFFE(2)-FFFS(2),FFFE(3)-FFFS(3) + 455 IF(JFT.LT.0)GO TO 475 + 460 CALL FFRT(AC,WS,FFTE,FFTS) + IF(JSS.NE.1)GO TO 470 + WRITE(67,*)FFTE(1),FFTS(1),FFTE(1)-FFTS(1) + WRITE(67,*)FFTE(2),FFTS(2),FFTE(2)-FFTS(2) + WRITE(67,*)FFTE(3),FFTS(3),FFTE(3)-FFTS(3) + GO TO 475 + 470 WRITE(ITT)(FFTE(I)-FFTS(I),I=1,3) + IF(JTW.EQ.1)WRITE(67,6666) + 1IX,IY,IZ,FFTE(1)-FFTS(1),FFTE(2)-FFTS(2),FFTE(3)-FFTS(3) + 475 CONTINUE + IF(JSS.EQ.1)GO TO 500 + IF(JFT.LE.0)CLOSE(ITF) + IF(JFT.GE.0)CLOSE(ITT) + 500 CONTINUE + CLOSE(IT) + STOP + 6666 FORMAT(I5,I4,I4,1PD15.4,1PD12.4,1PD12.4) + END + SUBROUTINE SAMPOA(AC,TMS,WS) + IMPLICIT REAL*8(A-H,O-Z) + COMMON/CIL/LE,NLEM,NLEMT,MXMPO,MXIM +CCC DIMENSION AC(NLEMT),TMS(LE,3),WS(NLEMT) + DIMENSION AC(880),TMS(20,3),WS(880) + COMPLEX*16 AC,TMS,WS + DIMENSION TM(2,2) + COMPLEX*16 TM + I=0 + DO 20 L=1,LE + TM(1,1)=TMS(L,1) + TM(1,2)=TMS(L,2) + TM(2,2)=TMS(L,3) + TM(2,1)=TM(1,2) + LTPO=L+L+1 + DO 20 IM=1,LTPO + I=I+1 + IE=I+NLEM + AC(I)=TM(1,1)*WS(I)+TM(1,2)*WS(IE) + AC(IE)=TM(2,1)*WS(I)+TM(2,2)*WS(IE) + 20 CONTINUE + RETURN + END + SUBROUTINE SAMP(AC,TMSM,TMSE,WS) + IMPLICIT REAL*8(A-H,O-Z) + COMMON/CIL/LE,NLEM,NLEMT,MXMPO,MXIM +CCC DIMENSION AC(NLEMT),TMSM(LE),TMSE(LE),WS(NLEMT) + DIMENSION AC(880),TMSM(20),TMSE(20),WS(880) + COMPLEX*16 AC,TMSM,TMSE,WS + I=0 + DO 20 L=1,LE + LTPO=L+L+1 + DO 20 IM=1,LTPO + I=I+1 + IE=I+NLEM + AC(I)=TMSM(L)*WS(I) + AC(IE)=TMSE(L)*WS(IE) + 20 CONTINUE + RETURN + END + SUBROUTINE CAMP(AC,AM0M,WS) + IMPLICIT REAL*8(A-H,O-Z) + COMMON/CIL/LE,NLEM,NLEMT,MXMPO,MXIM + COMMON/CCR/CC0,UIM,SQ2ITI,COF,CIMU,SQ2I + COMPLEX*16 CC0,UIM,SQ2ITI +CCC DIMENSION AC(NLEMT),AM0M(NLEMT,NLEMT),WS(NLEMT) + DIMENSION AC(880),AM0M(880,880),WS(880) + COMPLEX*16 AC,AM0M,WS + DO 20 J=1,NLEMT + AC(J)=CC0 + DO 20 I=1,NLEMT + AC(J)=AC(J)+AM0M(J,I)*WS(I) + 20 CONTINUE + RETURN + END + SUBROUTINE CZAMP(AC,AMD,INDAM,WS) + IMPLICIT REAL*8(A-H,O-Z) + COMMON/CIL/LE,NLEM,NLEMT,MXMPO,MXIM + COMMON/CCR/CC0,UIM,SQ2ITI,COF,CIMU,SQ2I + COMPLEX*16 CC0,UIM,SQ2ITI +CCC DIMENSION AC(NLEMT),AMD(NVAM,4),INDAM(LE,LE),WS(NLEMT) + DIMENSION AC(880),AMD(6140,4),INDAM(20,20),WS(880) + COMPLEX*16 AC,AMD,WS + COMPLEX*16 SUMM,SUME + DO 10 I=1,NLEMT + 10 AC(I)=CC0 + DO 20 IM=1,MXIM + M=IM-MXMPO + LMN=MAX0(IABS(M),1) + DO 20 L=LMN,LE + I=M+L*(L+1) + IE=I+NLEM + SUMM=CC0 + SUME=CC0 + DO 15 LS=LMN,LE + IS=M+LS*(LS+1) + ISE=IS+NLEM + NUM=INDAM(L,LS)+M + SUMM=SUMM+AMD(NUM,1)*WS(IS)+AMD(NUM,2)*WS(ISE) + SUME=SUME+AMD(NUM,3)*WS(IS)+AMD(NUM,4)*WS(ISE) + 15 CONTINUE + AC(I)=SUMM + AC(IE)=SUME + 20 CONTINUE + RETURN + END + SUBROUTINE FFRF(ZPV,AC,WS,FFFE,FFFS) + IMPLICIT REAL*8(A-H,O-Z) + COMMON/CIL/LE,NLEM,NLEMT,MXMPO,MXIM + COMMON/CCR/CC0,UIM,SQ2ITI,COF,CIMU,SQ2I + COMPLEX*16 CC0,UIM,SQ2ITI +CCC DIMENSION ZPV(LE,3,2,2),AC(NLEMT),WS(NLEMT) + DIMENSION ZPV(20,3,2,2),AC(880),WS(880),FFFE(3),FFFS(3) + COMPLEX*16 AC,WS + DIMENSION GAP(3) + COMPLEX*16 GAP + COMPLEX*16 UIMMP,SUMM,SUME,SUEM,SUEE + DO 50 IMU=1,3 + MU=IMU-2 + GAP(IMU)=CC0 + DO 40 L=1,LE + LPO=L+1 + LTPO=LPO+L + IMM=L*LPO + DO 40 ILMP=1,3 + IF(((L.EQ.1).AND.(ILMP.EQ.1)).OR.((L.EQ.LE).AND.(ILMP.EQ.3))) + 1GO TO 40 + LMPML=ILMP-2 + LMP=L+LMPML + UIMMP=-LMPML*UIM + IMPMMMP=LMP*(LMP+1) + DO 30 IM=1,LTPO + M=IM-LPO + MMP=M-MU + IF(IABS(MMP).GT.LMP)GO TO 30 + I=IMM+M + IE=I+NLEM + IMP=IMPMMMP+MMP + IMPE=IMP+NLEM + CGC=CG1(LMPML,MU,L,M) + SUMM=DCONJG(WS(I))*AC(IMP) + SUME=DCONJG(WS(I))*AC(IMPE) + SUEM=DCONJG(WS(IE))*AC(IMP) + SUEE=DCONJG(WS(IE))*AC(IMPE) + IF(LMPML.EQ.0)GO TO 25 + SUMM=SUMM*UIMMP + SUME=SUME*UIMMP + SUEM=SUEM*UIMMP + SUEE=SUEE*UIMMP + 25 GAP(IMU)=GAP(IMU)+CGC* + 1(SUMM*ZPV(L,ILMP,1,1)+SUME*ZPV(L,ILMP,1,2)+ + 2SUEM*ZPV(L,ILMP,2,1)+SUEE*ZPV(L,ILMP,2,2)) + 30 CONTINUE + 40 CONTINUE + 50 CONTINUE + SUME=-GAP(1)*CIMU + SUEE=-GAP(2)*COF + SUEM=-GAP(3)*CIMU + FFFE(1)=DREAL(SUME-SUEM) + FFFE(2)=DREAL((SUME+SUEM)*UIM) + FFFE(3)=DREAL(SUEE) + DO 90 IMU=1,3 + MU=IMU-2 + GAP(IMU)=CC0 + DO 80 L=1,LE + LPO=L+1 + LTPO=LPO+L + IMM=L*LPO + DO 80 ILMP=1,3 + IF(((L.EQ.1).AND.(ILMP.EQ.1)).OR.((L.EQ.LE).AND.(ILMP.EQ.3))) + 1GO TO 80 + LMPML=ILMP-2 + LMP=L+LMPML + UIMMP=-LMPML*UIM + IMPMMMP=LMP*(LMP+1) + DO 70 IM=1,LTPO + M=IM-LPO + MMP=M-MU + IF(IABS(MMP).GT.LMP)GO TO 70 + I=IMM+M + IE=I+NLEM + IMP=IMPMMMP+MMP + IMPE=IMP+NLEM + CGC=CG1(LMPML,MU,L,M) + SUMM=DCONJG(AC(I))*AC(IMP) + SUME=DCONJG(AC(I))*AC(IMPE) + SUEM=DCONJG(AC(IE))*AC(IMP) + SUEE=DCONJG(AC(IE))*AC(IMPE) + IF(LMPML.EQ.0)GO TO 65 + SUMM=SUMM*UIMMP + SUME=SUME*UIMMP + SUEM=SUEM*UIMMP + SUEE=SUEE*UIMMP + 65 GAP(IMU)=GAP(IMU)+CGC* + 1(SUMM*ZPV(L,ILMP,1,1)+SUME*ZPV(L,ILMP,1,2)+ + 2SUEM*ZPV(L,ILMP,2,1)+SUEE*ZPV(L,ILMP,2,2)) + 70 CONTINUE + 80 CONTINUE + 90 CONTINUE + SUME=GAP(1)*CIMU + SUEE=GAP(2)*COF + SUEM=GAP(3)*CIMU + FFFS(1)=DREAL(SUME-SUEM) + FFFS(2)=DREAL((SUME+SUEM)*UIM) + FFFS(3)=DREAL(SUEE) + RETURN + END + SUBROUTINE THDPS(LM,ZPV) + IMPLICIT REAL*8(A-H,O-Z) +CCC DIMENSION ZPV(LM,3,2,2) + DIMENSION ZPV(20,3,2,2) + C0=0.0D0 + DO 10 L=1,LM + DO 10 ILMP=1,3 + ZPV(L,ILMP,1,1)=C0 + ZPV(L,ILMP,1,2)=C0 + ZPV(L,ILMP,2,1)=C0 + ZPV(L,ILMP,2,2)=C0 + 10 CONTINUE + DO 15 L=1,LM + XD=L*(L+1) + ZP=-1.0D0/DSQRT(XD) + ZPV(L,2,1,2)=ZP + ZPV(L,2,2,1)=ZP + 15 CONTINUE + IF(LM.EQ.1)GO TO 30 + DO 20 L=2,LM + XN=(L-1)*(L+1) + XD=L*(L+L+1) + ZP=DSQRT(XN/XD) + ZPV(L,1,1,1)=ZP + ZPV(L,1,2,2)=ZP + 20 CONTINUE + LMMO=LM-1 + DO 25 L=1,LMMO + XN=L*(L+2) + XD=(L+1)*(L+L+1) + ZP=-DSQRT(XN/XD) + ZPV(L,3,1,1)=ZP + ZPV(L,3,2,2)=ZP + 25 CONTINUE + 30 CONTINUE + RETURN + END + REAL*8 FUNCTION CG1(LMPML,MU,L,M) +CCC CG1(LMPML,MU,L,M)=CLGO(1,LMP,L;MU,M-MU,M) + REAL*8 XD,XN + IF(LMPML)30,5,60 + 5 IF((M.NE.0).OR.(MU.NE.0))GO TO 10 + CG1=0.0D0 + RETURN + 10 IF(MU.EQ.0)GO TO 20 + XD=(L+1)*L*2 + IF(MU.GT.0)GO TO 15 + XN=(L-M)*(L+M+1) + CG1=-DSQRT(XN/XD) + RETURN + 15 XN=(L+M)*(L-M+1) + CG1=DSQRT(XN/XD) + RETURN + 20 XD=(L+1)*L + XN=-M + CG1=XN/DSQRT(XD) + RETURN + 30 XD=(L*2-1)*L*2 + IF(MU)35,40,45 + 35 XN=(L-1-M)*(L-M) + CG1=DSQRT(XN/XD) + RETURN + 40 XN=(L-M)*(L+M)*2 + CG1=DSQRT(XN/XD) + RETURN + 45 XN=(L-1+M)*(L+M) + CG1=DSQRT(XN/XD) + RETURN + 60 XD=(L*2+3)*(L+1)*2 + IF(MU)65,70,75 + 65 XN=(L+1+M)*(L+2+M) + CG1=DSQRT(XN/XD) + RETURN + 70 XN=(L+1-M)*(L+1+M)*2 + CG1=-DSQRT(XN/XD) + RETURN + 75 XN=(L+1-M)*(L+2-M) + CG1=DSQRT(XN/XD) + RETURN + END + SUBROUTINE FFRT(AC,WS,FFTE,FFTS) + IMPLICIT REAL*8(A-H,O-Z) + COMMON/CIL/LE,NLEM,NLEMT,MXMPO,MXIM + COMMON/CCR/CC0,UIM,SQ2ITI,COF,CIMU,SQ2I + COMPLEX*16 CC0,UIM,SQ2ITI +CCC DIMENSION AC(NLEMT),WS(NLEMT) + DIMENSION AC(880),WS(880),FFTE(3),FFTS(3) + COMPLEX*16 AC,WS + DIMENSION CTQCE(3),CTQCS(3) + COMPLEX*16 CTQCE,CTQCS + COMPLEX*16 ACW,ACA + CTQCE(1)=CC0 + CTQCE(2)=CC0 + CTQCE(3)=CC0 + CTQCS(1)=CC0 + CTQCS(2)=CC0 + CTQCS(3)=CC0 + DO 60 L=1,LE + LPO=L+1 + IL=L*LPO + LTPO=L+LPO + DO 60 IM=1,LTPO + M=IM-LPO + I=M+IL + IE=I+NLEM + MMMU=M+1 + MMMMU=IABS(MMMU) + IF(MMMMU.GT.L)GO TO 30 + IMMU=MMMU+IL + IMMUE=IMMU+NLEM + RMU=(L+MMMU)*(L-M) + RMU=-DSQRT(RMU)*SQ2I + ACW=DCONJG(AC(I))*WS(IMMU)+DCONJG(AC(IE))*WS(IMMUE) + ACA=DCONJG(AC(I))*AC(IMMU)+DCONJG(AC(IE))*AC(IMMUE) + CTQCE(1)=CTQCE(1)+ACW*RMU + CTQCS(1)=CTQCS(1)+ACA*RMU + 30 RMU=-M + ACW=DCONJG(AC(I))*WS(I)+DCONJG(AC(IE))*WS(IE) + ACA=DCONJG(AC(I))*AC(I)+DCONJG(AC(IE))*AC(IE) + CTQCE(2)=CTQCE(2)+ACW*RMU + CTQCS(2)=CTQCS(2)+ACA*RMU + MMMU=M-1 + MMMMU=IABS(MMMU) + IF(MMMMU.GT.L)GO TO 50 + IMMU=MMMU+IL + IMMUE=IMMU+NLEM + RMU=(L-MMMU)*(L+M) + RMU=DSQRT(RMU)*SQ2I + ACW=DCONJG(AC(I))*WS(IMMU)+DCONJG(AC(IE))*WS(IMMUE) + ACA=DCONJG(AC(I))*AC(IMMU)+DCONJG(AC(IE))*AC(IMMUE) + CTQCE(3)=CTQCE(3)+ACW*RMU + CTQCS(3)=CTQCS(3)+ACA*RMU + 50 CONTINUE + 60 CONTINUE + FFTE(1)=DREAL(CTQCE(1)-CTQCE(3))*SQ2I + FFTE(2)=DREAL((CTQCE(1)+CTQCE(3))*SQ2ITI) + FFTE(3)=DREAL(CTQCE(2)) + FFTS(1)=-DREAL(CTQCS(1)-CTQCS(3))*SQ2I + FFTS(2)=-DREAL((CTQCS(1)+CTQCS(3))*SQ2ITI) + FFTS(3)=-DREAL(CTQCS(2)) + RETURN + END +CCC \ No newline at end of file -- GitLab