Upload New File

trapping/lffft.f

+      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