Loading src/cluster/clu.f +2 −29 Original line number Diff line number Diff line Loading @@ -330,12 +330,8 @@ C 1", arg =",ARG GO TO 490 132 CONTINUE CALL CMS(AM) CALL SUMMAT(AM,CCSAM) PRINT *,"DEBUG: after CMS CCSAM =",CCSAM NDIT=2*NSPH*NLIM CALL LUCIN(AM,MXNDM,NDIT,JER) CALL SUMMAT(AM,CCSAM) PRINT *,"DEBUG: after LUCIN CCSAM =",CCSAM IF(JER.EQ.1)GO TO 495 CALL ZTM(AM) Loading @@ -358,7 +354,6 @@ C 158 WRITE(IW,6061) 160 CS0=0.25D0*VK*VK*VK/DACOS(C0) CALL SCR0(VK,EXRI) PRINT *,"DEBUG: after SCR0 TFSAS =",TFSAS SQK=VK*VK*EXDC CALL APS(ZPV,LI,NSPH,IOG,RMI,REI,SQK,GAPS) CALL RABAS(INPOL,LI,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS) Loading Loading @@ -639,6 +634,7 @@ C 482 IF(ISAM.LE.1)THSL=THSL+THSSTP 484 PH=PH+PHSTP 486 TH=TH+THSTP PRINT *,"INFO: done scale iteration ",JXI 488 CONTINUE GO TO 500 490 WRITE(IW,6550)JER,LCALC,ARG Loading @@ -650,6 +646,7 @@ C CLOSE(IW) CLOSE(IT) CLOSE(ITIN) PRINT *,"INFO: output written to OCLU." STOP END SUBROUTINE RABAS(INPOL,LI,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS) Loading Loading @@ -3504,28 +3501,4 @@ CCC DIMENSION AV(NDDMST*NDDMST) 3 CDTP=CDTP+AV((J-1)*ISTEP+I)*AV(J+IR) RETURN END SUBROUTINE SUMMAT(AM,RSLT) DIMENSION AM(921600) COMPLEX*16 AM,CC0,RSLT CC0=(0.0D0,0.0D0) RSLT=CC0 DO 4 I=1,921600 4 RSLT=RSLT+AM(I) RETURN END SUBROUTINE LOGMAT(AM) DIMENSION AM(960,960) COMPLEX*16 AM IL=36 OPEN(IL,FILE="matrix.log",STATUS="UNKNOWN") DO 5 I=1,960 DO 6 J=1,960 WRITE(IL,9901)I,J,DREAL(AM(I,J)) 6 WRITE(IL,9902)I,J,DIMAG(AM(I,J)) 5 CONTINUE CLOSE(IL) 9901 FORMAT("R:",1I3,",",1I3,",",1D15.7) 9902 FORMAT("I:",1I3,",",1I3,",",1D15.7) RETURN END CCC src/sphere/sph.f +2 −61 Original line number Diff line number Diff line Loading @@ -270,7 +270,6 @@ C CS0=0.25D0*VK*VK*VK/PIGH CALL SSCR0(TFSAS,NSPH,LM,VK,EXRI) PRINT *,"DEBUG: TFSAS =", TFSAS SQK=VK*VK*EXDC CALL APS(ZPV,LM,NSPH,IOG,RMI,REI,SQK,GAPS) CALL RABAS(INPOL,LM,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS) Loading Loading @@ -397,6 +396,7 @@ C 482 IF(ISAM.LE.1)THSL=THSL+THSSTP 484 PH=PH+PHSTP 486 TH=TH+THSTP PRINT *,"INFO: done scale iteration ",JXI 488 CONTINUE GO TO 500 490 WRITE(IW,6550)JER,LCALC,ARG Loading @@ -406,6 +406,7 @@ C CLOSE(IW) CLOSE(IT) CLOSE(ITIN) PRINT *,"INFO: output written to OSPH." STOP END SUBROUTINE RABAS(INPOL,LI,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS) Loading Loading @@ -456,19 +457,11 @@ CCC 1TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH) TQSPE(2,I)=TQSPE(2,I)*PIG2 TQSPS(1,I)=TQSPS(1,I)*PIG2 TQSPS(2,I)=TQSPS(2,I)*PIG2 PRINT *,"DEBUG: TQSPE(1,",I,") =",TQSPE(1,I) PRINT *,"DEBUG: TQSPE(2,",I,") =",TQSPE(2,I) PRINT *,"DEBUG: TQSPS(1,",I,") =",TQSPS(1,I) PRINT *,"DEBUG: TQSPS(2,",I,") =",TQSPS(2,I) GO TO 80 75 TQSE(1,I)=TQSE(1,I)*PIG2 TQSE(2,I)=TQSE(2,I)*PIG2 TQSS(1,I)=TQSS(1,I)*PIG2 TQSS(2,I)=TQSS(2,I)*PIG2 PRINT *,"DEBUG: TQSE(1,",I,") =",TQSE(1,I) PRINT *,"DEBUG: TQSE(2,",I,") =",TQSE(2,I) PRINT *,"DEBUG: TQSS(1,",I,") =",TQSS(1,I) PRINT *,"DEBUG: TQSS(2,",I,") =",TQSS(2,I) 80 CONTINUE RETURN END Loading Loading @@ -584,7 +577,6 @@ CCC 1TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH) PIGFSQ=6.4D+1*DACOS(C0)**2 CFSQ=4.0D0/(PIGFSQ*CCS*CCS) NLMM=LM*(LM+2) C PRINT *,"DEBUG: in SSCR2 W(1,1) =",W(1,1) DO 14 I=1,NSPH IOGI=IOG(I) IF(IOGI.LT.I)GO TO 14 Loading @@ -596,30 +588,14 @@ C PRINT *,"DEBUG: in SSCR2 W(1,1) =",W(1,1) DO 10 L=1,LM RM=1.0D0/RMI(L,I) RE=1.0D0/REI(L,I) C PRINT *,"DEBUG: RM =",RM C PRINT *,"DEBUG: RE =",RE LTPO=L+L+1 DO 10 IM=1,LTPO K=K+1 KE=K+NLMM C PRINT *,"DEBUG: W(",K,",3) =",W(K,3) C PRINT *,"DEBUG: W(",K,",1) =",W(K,1) C PRINT *,"DEBUG: W(",KE,",3) =",W(KE,3) C PRINT *,"DEBUG: W(",KE,",1) =",W(KE,1) S11=S11-W(K,3)*W(K,1)*RM-W(KE,3)*W(KE,1)*RE C PRINT *,"DEBUG: S11 =",S11 C PRINT *,"DEBUG: W(",K,",4) =",W(K,4) C PRINT *,"DEBUG: W(",KE,",4) =",W(KE,4) S21=S21-W(K,4)*W(K,1)*RM-W(KE,4)*W(KE,1)*RE C PRINT *,"DEBUG: W(",K,",2) =",W(K,2) C PRINT *,"DEBUG: W(",KE,",2) =",W(KE,2) S12=S12-W(K,3)*W(K,2)*RM-W(KE,3)*W(KE,2)*RE 10 S22=S22-W(K,4)*W(K,2)*RM-W(KE,4)*W(KE,2)*RE C PRINT *,"DEBUG: CSAM =",CSAM C PRINT *,"DEBUG: S11 =",S11 C PRINT *,"DEBUG: S21 =",S21 C PRINT *,"DEBUG: S12 =",S12 C PRINT *,"DEBUG: S22 =",S22 SAS(I,1,1)=S11*CSAM SAS(I,2,1)=S21*CSAM SAS(I,1,2)=S12*CSAM Loading Loading @@ -670,7 +646,6 @@ CCC 1RMI(LI,NSPH),REI(LI,NSPH),GAPS(NSPH) 1CG1(LMPML,0,L,1)*(SUMM+SUME+SUEM+SUEE))*COFL 30 CONTINUE GAPS(I)=DREAL(SUM)*COFS PRINT *,"DEBUG: GAPS(",I,") =",GAPS(I) 40 CONTINUE RETURN END Loading Loading @@ -781,10 +756,6 @@ CCC CG1(LMPML,MU,L,M)=CLGO(1,LMP,L;MU,M-MU,M) SINT=DSIN(TH) COSP=DCOS(PH) SINP=DSIN(PH) C PRINT *,"DEBUG: cost =",COST C PRINT *,"DEBUG: sint =",SINT C PRINT *,"DEBUG: cosp =",COSP C PRINT *,"DEBUG: sinp =",SINP U(1)=COSP*SINT U(2)=SINP*SINT U(3)=COST Loading Loading @@ -828,8 +799,6 @@ CCC NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE DO 60 N=1,NSPH 60 ARG(N)=-ARG(N) 65 CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM) C PRINT *,"DEBUG: in WMAMP and calling PWMA with lm =",LM, C 1"and iis =",IIS CALL PWMA(UP,UN,YLM,INPOL,LM,IIS) RETURN END Loading Loading @@ -928,11 +897,9 @@ CCC NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE 55 ARGS(N)=-COSTS*RZZ(N) 60 CONTINUE 75 CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM) C PRINT *,"DEBUG: in WMASP and calling PWMA" CALL PWMA(UP,UN,YLM,INPOL,LM,ISQ) IF(IBF.LT.0)RETURN CALL SPHAR(COSTS,SINTS,COSPS,SINPS,LM,YLM) C PRINT *,"DEBUG: in WMASP and calling PWMA with IBF =",IBF CALL PWMA(UPS,UNS,YLM,INPOL,LM,2) RETURN END Loading Loading @@ -963,7 +930,6 @@ CCC DIMENSION YLM(NLWM+2) CM2=.5D0*DCMPLX(UN(1),UN(2)) CP2=.5D0*DCMPLX(UN(1),-UN(2)) CZ2=UN(3) C PRINT *,"DEBUG: in PWMA YLM(2) =",YLM(2) DO 20 L=1,LW LF=L+1 LFTL=LF*L Loading @@ -980,35 +946,26 @@ C PRINT *,"DEBUG: in PWMA YLM(2) =",YLM(2) CM=DSQRT(X) CZ=M W(K,ISPO)=DCONJG(CP1*CP*YLM(K+2)+CM1*CM*YLM(K)+CZ1*CZ*YLM(K+1))*CL C IF(ISPO.EQ.1)PRINT *,"DEBUG: W(",K,",",ISPO,") =",W(K,ISPO) 20 W(K,ISPT)=DCONJG(CP2*CP*YLM(K+2)+CM2*CM*YLM(K)+CZ2*CZ*YLM(K+1))*CL C 20 PRINT *,"DEBUG: W(",K,",",ISPT,") =",W(K,ISPT) DO 30 K=1,NLWM I=K+NLWM W(I,ISPO)=UIM*W(K,ISPT) C PRINT *,"DEBUG: W(",I,",",ISPO,") =",W(I,ISPO) 30 W(I,ISPT)=-UIM*W(K,ISPO) C 30 PRINT *,"DEBUG: W(",I,",",ISPT,") =",W(I,ISPT) IF(INPOL.EQ.0)GO TO 42 DO 40 K=1,NLWM I=K+NLWM C1=(W(K,ISPO)+UIM*W(K,ISPT))*SQRTWI C2=(W(K,ISPO)-UIM*W(K,ISPT))*SQRTWI W(K,ISPO)=C2 C PRINT *,"DEBUG: W(",K,",",ISPO,") =",W(K,ISPO) W(I,ISPO)=-C2 C PRINT *,"DEBUG: W(",I,",",ISPO,") =",W(I,ISPO) W(K,ISPT)=C1 C PRINT *,"DEBUG: W(",K,",",ISPT,") =",W(K,ISPT) 40 W(I,ISPT)=C1 C 40 PRINT *,"DEBUG: W(",I,",",ISPT,") =",W(I,ISPT) 42 IF(ISQ.EQ.0)RETURN DO 50 I=1,2 IPT=I+2 IPIS=I+IS DO 50 K=1,NLWMT 50 W(K,IPT)=DCONJG(W(K,IPIS)) C 50 PRINT *,"DEBUG: W(",K,",",IPT,") =",W(K,IPT) RETURN END SUBROUTINE ORUNVE(U1,U2,U3,IORTH,TORTH) Loading Loading @@ -1101,9 +1058,7 @@ CCC 2RMF(LI),DRMF(LI),REF(LI),DREF(LI) CCNC=FBI(LPO)*DFB CCND=FB(LPO)*DFBI RMI(L,I)=1.0D0+UIM*(CCNA-CCNB)/(CCNC-CCND) C PRINT *,"DEBUG: gone 60, RMI(",L,",",I,") =",RMI(L,I) 60 REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND) C 60 PRINT *,"DEBUG: gone 60, REI(",L,",",I,") =",REI(L,I) RETURN 65 DO 80 L=1,LI LPO=L+1 Loading Loading @@ -1138,7 +1093,6 @@ cccccccccccccccccccccc 80 DREF(L)=DY2*SZ+Y2 CRI=(1.0D0,0.0D0) IF(MOD(NSH,2).NE.0)CRI=DC0(IC)/EXDC C PRINT *,"DEBUG: going 90, AREX =",AREX DO 90 L=1,LI LPO=L+1 LTPO=LPO+L Loading @@ -1150,13 +1104,11 @@ C PRINT *,"DEBUG: going 90, AREX =",AREX CCNC=RMF(L)*DFB CCND=DRMF(L)*FB(LPO)*SZ*LTPO RMI(L,I)=1.0D0+UIM*(CCNA-CCNB)/(CCNC-CCND) C PRINT *,"DEBUG: gone 90, RMI(",L,",",I,") =",RMI(L,I) CCNA=REF(L)*DFN CCNB=DREF(L)*FN(LPO)*SZ*LTPO CCNC=REF(L)*DFB CCND=DREF(L)*FB(LPO)*SZ*LTPO 90 REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND) C 90 PRINT *,"DEBUG: gone 90, REI(",L,",",I,") =",REI(L,I) RETURN END SUBROUTINE RKT(NPNTMO,STEP,X,LPO,Y1,Y2,DY1,DY2) Loading Loading @@ -1498,18 +1450,13 @@ CCC LL=LM PI4=DACOS(0.0D0)*8.0D0 PI4IRS=1.0D0/DSQRT(PI4) X=COSRTH C PRINT *,"DEBUG: X =",X Y=DABS(SINRTH) C PRINT *,"DEBUG: Y =",Y CLLMO=3.0D0 CLL=1.5D0 YTOL=Y PLEGN(1)=1.0D0 C PRINT *,"DEBUG: PLEGN( 1 ) =",PLEGN(1) PLEGN(2)=X*DSQRT(CLLMO) C PRINT *,"DEBUG: PLEGN( 2 ) =",PLEGN(2) PLEGN(3)=YTOL*DSQRT(CLL) C PRINT *,"DEBUG: PLEGN( 3 ) =",PLEGN(3) SINRMP(1)=SINRPH COSRMP(1)=COSRPH IF(LL.LT.2)GO TO 30 Loading @@ -1529,17 +1476,14 @@ C PRINT *,"DEBUG: PLEGN( 3 ) =",PLEGN(3) CDM=(LTMO-2)*LS 10 PLEGN(MPOPK)=PLEGN(MPOPK-L)*X*DSQRT(CN/CD)- 1PLEGN(MPOPK-LTMO)*DSQRT(CNM/CDM) C10 PRINT *,"DEBUG: PLEGN(",MPOPK,") =",PLEGN(MPOPK) LPK=L+K CLTPO=LTPO PLEGN(LPK)=PLEGN(K)*X*DSQRT(CLTPO) C PRINT *,"DEBUG: PLEGN(",LPK,") =",PLEGN(LPK) K=LPK+1 CLT=LTPO-1 CLL=CLL*(CLTPO/CLT) YTOL=YTOL*Y PLEGN(K)=YTOL*DSQRT(CLL) C PRINT *,"DEBUG: PLEGN(",K,") =",PLEGN(K) SINRMP(L)=SINRPH*COSRMP(LMO)+COSRPH*SINRMP(LMO) 20 COSRMP(L)=COSRPH*COSRMP(LMO)-SINRPH*SINRMP(LMO) 30 L=0 Loading @@ -1548,17 +1492,14 @@ C PRINT *,"DEBUG: PLEGN(",K,") =",PLEGN(K) L0Y=K+1 L0P=K/2+1 YLM(L0Y)=PI4IRS*PLEGN(L0P) C PRINT *, "DEBUG: YLM(",L0Y,") =",YLM(L0Y) 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) C PRINT *, "DEBUG: YLM(",LMY,") =",YLM(LMY) LMY=L0Y-M YLM(LMY)=SAVE*DCMPLX(COSRMP(M),-SINRMP(M)) C PRINT *, "DEBUG: YLM(",LMY,") =",YLM(LMY) 45 IF(M.GE.L)GO TO 47 M=M+1 GO TO 44 Loading src/sphere/sphere.cpp +1 −1 File changed.Contains only whitespace changes. Show changes Loading
src/cluster/clu.f +2 −29 Original line number Diff line number Diff line Loading @@ -330,12 +330,8 @@ C 1", arg =",ARG GO TO 490 132 CONTINUE CALL CMS(AM) CALL SUMMAT(AM,CCSAM) PRINT *,"DEBUG: after CMS CCSAM =",CCSAM NDIT=2*NSPH*NLIM CALL LUCIN(AM,MXNDM,NDIT,JER) CALL SUMMAT(AM,CCSAM) PRINT *,"DEBUG: after LUCIN CCSAM =",CCSAM IF(JER.EQ.1)GO TO 495 CALL ZTM(AM) Loading @@ -358,7 +354,6 @@ C 158 WRITE(IW,6061) 160 CS0=0.25D0*VK*VK*VK/DACOS(C0) CALL SCR0(VK,EXRI) PRINT *,"DEBUG: after SCR0 TFSAS =",TFSAS SQK=VK*VK*EXDC CALL APS(ZPV,LI,NSPH,IOG,RMI,REI,SQK,GAPS) CALL RABAS(INPOL,LI,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS) Loading Loading @@ -639,6 +634,7 @@ C 482 IF(ISAM.LE.1)THSL=THSL+THSSTP 484 PH=PH+PHSTP 486 TH=TH+THSTP PRINT *,"INFO: done scale iteration ",JXI 488 CONTINUE GO TO 500 490 WRITE(IW,6550)JER,LCALC,ARG Loading @@ -650,6 +646,7 @@ C CLOSE(IW) CLOSE(IT) CLOSE(ITIN) PRINT *,"INFO: output written to OCLU." STOP END SUBROUTINE RABAS(INPOL,LI,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS) Loading Loading @@ -3504,28 +3501,4 @@ CCC DIMENSION AV(NDDMST*NDDMST) 3 CDTP=CDTP+AV((J-1)*ISTEP+I)*AV(J+IR) RETURN END SUBROUTINE SUMMAT(AM,RSLT) DIMENSION AM(921600) COMPLEX*16 AM,CC0,RSLT CC0=(0.0D0,0.0D0) RSLT=CC0 DO 4 I=1,921600 4 RSLT=RSLT+AM(I) RETURN END SUBROUTINE LOGMAT(AM) DIMENSION AM(960,960) COMPLEX*16 AM IL=36 OPEN(IL,FILE="matrix.log",STATUS="UNKNOWN") DO 5 I=1,960 DO 6 J=1,960 WRITE(IL,9901)I,J,DREAL(AM(I,J)) 6 WRITE(IL,9902)I,J,DIMAG(AM(I,J)) 5 CONTINUE CLOSE(IL) 9901 FORMAT("R:",1I3,",",1I3,",",1D15.7) 9902 FORMAT("I:",1I3,",",1I3,",",1D15.7) RETURN END CCC
src/sphere/sph.f +2 −61 Original line number Diff line number Diff line Loading @@ -270,7 +270,6 @@ C CS0=0.25D0*VK*VK*VK/PIGH CALL SSCR0(TFSAS,NSPH,LM,VK,EXRI) PRINT *,"DEBUG: TFSAS =", TFSAS SQK=VK*VK*EXDC CALL APS(ZPV,LM,NSPH,IOG,RMI,REI,SQK,GAPS) CALL RABAS(INPOL,LM,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS) Loading Loading @@ -397,6 +396,7 @@ C 482 IF(ISAM.LE.1)THSL=THSL+THSSTP 484 PH=PH+PHSTP 486 TH=TH+THSTP PRINT *,"INFO: done scale iteration ",JXI 488 CONTINUE GO TO 500 490 WRITE(IW,6550)JER,LCALC,ARG Loading @@ -406,6 +406,7 @@ C CLOSE(IW) CLOSE(IT) CLOSE(ITIN) PRINT *,"INFO: output written to OSPH." STOP END SUBROUTINE RABAS(INPOL,LI,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS) Loading Loading @@ -456,19 +457,11 @@ CCC 1TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH) TQSPE(2,I)=TQSPE(2,I)*PIG2 TQSPS(1,I)=TQSPS(1,I)*PIG2 TQSPS(2,I)=TQSPS(2,I)*PIG2 PRINT *,"DEBUG: TQSPE(1,",I,") =",TQSPE(1,I) PRINT *,"DEBUG: TQSPE(2,",I,") =",TQSPE(2,I) PRINT *,"DEBUG: TQSPS(1,",I,") =",TQSPS(1,I) PRINT *,"DEBUG: TQSPS(2,",I,") =",TQSPS(2,I) GO TO 80 75 TQSE(1,I)=TQSE(1,I)*PIG2 TQSE(2,I)=TQSE(2,I)*PIG2 TQSS(1,I)=TQSS(1,I)*PIG2 TQSS(2,I)=TQSS(2,I)*PIG2 PRINT *,"DEBUG: TQSE(1,",I,") =",TQSE(1,I) PRINT *,"DEBUG: TQSE(2,",I,") =",TQSE(2,I) PRINT *,"DEBUG: TQSS(1,",I,") =",TQSS(1,I) PRINT *,"DEBUG: TQSS(2,",I,") =",TQSS(2,I) 80 CONTINUE RETURN END Loading Loading @@ -584,7 +577,6 @@ CCC 1TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH) PIGFSQ=6.4D+1*DACOS(C0)**2 CFSQ=4.0D0/(PIGFSQ*CCS*CCS) NLMM=LM*(LM+2) C PRINT *,"DEBUG: in SSCR2 W(1,1) =",W(1,1) DO 14 I=1,NSPH IOGI=IOG(I) IF(IOGI.LT.I)GO TO 14 Loading @@ -596,30 +588,14 @@ C PRINT *,"DEBUG: in SSCR2 W(1,1) =",W(1,1) DO 10 L=1,LM RM=1.0D0/RMI(L,I) RE=1.0D0/REI(L,I) C PRINT *,"DEBUG: RM =",RM C PRINT *,"DEBUG: RE =",RE LTPO=L+L+1 DO 10 IM=1,LTPO K=K+1 KE=K+NLMM C PRINT *,"DEBUG: W(",K,",3) =",W(K,3) C PRINT *,"DEBUG: W(",K,",1) =",W(K,1) C PRINT *,"DEBUG: W(",KE,",3) =",W(KE,3) C PRINT *,"DEBUG: W(",KE,",1) =",W(KE,1) S11=S11-W(K,3)*W(K,1)*RM-W(KE,3)*W(KE,1)*RE C PRINT *,"DEBUG: S11 =",S11 C PRINT *,"DEBUG: W(",K,",4) =",W(K,4) C PRINT *,"DEBUG: W(",KE,",4) =",W(KE,4) S21=S21-W(K,4)*W(K,1)*RM-W(KE,4)*W(KE,1)*RE C PRINT *,"DEBUG: W(",K,",2) =",W(K,2) C PRINT *,"DEBUG: W(",KE,",2) =",W(KE,2) S12=S12-W(K,3)*W(K,2)*RM-W(KE,3)*W(KE,2)*RE 10 S22=S22-W(K,4)*W(K,2)*RM-W(KE,4)*W(KE,2)*RE C PRINT *,"DEBUG: CSAM =",CSAM C PRINT *,"DEBUG: S11 =",S11 C PRINT *,"DEBUG: S21 =",S21 C PRINT *,"DEBUG: S12 =",S12 C PRINT *,"DEBUG: S22 =",S22 SAS(I,1,1)=S11*CSAM SAS(I,2,1)=S21*CSAM SAS(I,1,2)=S12*CSAM Loading Loading @@ -670,7 +646,6 @@ CCC 1RMI(LI,NSPH),REI(LI,NSPH),GAPS(NSPH) 1CG1(LMPML,0,L,1)*(SUMM+SUME+SUEM+SUEE))*COFL 30 CONTINUE GAPS(I)=DREAL(SUM)*COFS PRINT *,"DEBUG: GAPS(",I,") =",GAPS(I) 40 CONTINUE RETURN END Loading Loading @@ -781,10 +756,6 @@ CCC CG1(LMPML,MU,L,M)=CLGO(1,LMP,L;MU,M-MU,M) SINT=DSIN(TH) COSP=DCOS(PH) SINP=DSIN(PH) C PRINT *,"DEBUG: cost =",COST C PRINT *,"DEBUG: sint =",SINT C PRINT *,"DEBUG: cosp =",COSP C PRINT *,"DEBUG: sinp =",SINP U(1)=COSP*SINT U(2)=SINP*SINT U(3)=COST Loading Loading @@ -828,8 +799,6 @@ CCC NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE DO 60 N=1,NSPH 60 ARG(N)=-ARG(N) 65 CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM) C PRINT *,"DEBUG: in WMAMP and calling PWMA with lm =",LM, C 1"and iis =",IIS CALL PWMA(UP,UN,YLM,INPOL,LM,IIS) RETURN END Loading Loading @@ -928,11 +897,9 @@ CCC NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE 55 ARGS(N)=-COSTS*RZZ(N) 60 CONTINUE 75 CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM) C PRINT *,"DEBUG: in WMASP and calling PWMA" CALL PWMA(UP,UN,YLM,INPOL,LM,ISQ) IF(IBF.LT.0)RETURN CALL SPHAR(COSTS,SINTS,COSPS,SINPS,LM,YLM) C PRINT *,"DEBUG: in WMASP and calling PWMA with IBF =",IBF CALL PWMA(UPS,UNS,YLM,INPOL,LM,2) RETURN END Loading Loading @@ -963,7 +930,6 @@ CCC DIMENSION YLM(NLWM+2) CM2=.5D0*DCMPLX(UN(1),UN(2)) CP2=.5D0*DCMPLX(UN(1),-UN(2)) CZ2=UN(3) C PRINT *,"DEBUG: in PWMA YLM(2) =",YLM(2) DO 20 L=1,LW LF=L+1 LFTL=LF*L Loading @@ -980,35 +946,26 @@ C PRINT *,"DEBUG: in PWMA YLM(2) =",YLM(2) CM=DSQRT(X) CZ=M W(K,ISPO)=DCONJG(CP1*CP*YLM(K+2)+CM1*CM*YLM(K)+CZ1*CZ*YLM(K+1))*CL C IF(ISPO.EQ.1)PRINT *,"DEBUG: W(",K,",",ISPO,") =",W(K,ISPO) 20 W(K,ISPT)=DCONJG(CP2*CP*YLM(K+2)+CM2*CM*YLM(K)+CZ2*CZ*YLM(K+1))*CL C 20 PRINT *,"DEBUG: W(",K,",",ISPT,") =",W(K,ISPT) DO 30 K=1,NLWM I=K+NLWM W(I,ISPO)=UIM*W(K,ISPT) C PRINT *,"DEBUG: W(",I,",",ISPO,") =",W(I,ISPO) 30 W(I,ISPT)=-UIM*W(K,ISPO) C 30 PRINT *,"DEBUG: W(",I,",",ISPT,") =",W(I,ISPT) IF(INPOL.EQ.0)GO TO 42 DO 40 K=1,NLWM I=K+NLWM C1=(W(K,ISPO)+UIM*W(K,ISPT))*SQRTWI C2=(W(K,ISPO)-UIM*W(K,ISPT))*SQRTWI W(K,ISPO)=C2 C PRINT *,"DEBUG: W(",K,",",ISPO,") =",W(K,ISPO) W(I,ISPO)=-C2 C PRINT *,"DEBUG: W(",I,",",ISPO,") =",W(I,ISPO) W(K,ISPT)=C1 C PRINT *,"DEBUG: W(",K,",",ISPT,") =",W(K,ISPT) 40 W(I,ISPT)=C1 C 40 PRINT *,"DEBUG: W(",I,",",ISPT,") =",W(I,ISPT) 42 IF(ISQ.EQ.0)RETURN DO 50 I=1,2 IPT=I+2 IPIS=I+IS DO 50 K=1,NLWMT 50 W(K,IPT)=DCONJG(W(K,IPIS)) C 50 PRINT *,"DEBUG: W(",K,",",IPT,") =",W(K,IPT) RETURN END SUBROUTINE ORUNVE(U1,U2,U3,IORTH,TORTH) Loading Loading @@ -1101,9 +1058,7 @@ CCC 2RMF(LI),DRMF(LI),REF(LI),DREF(LI) CCNC=FBI(LPO)*DFB CCND=FB(LPO)*DFBI RMI(L,I)=1.0D0+UIM*(CCNA-CCNB)/(CCNC-CCND) C PRINT *,"DEBUG: gone 60, RMI(",L,",",I,") =",RMI(L,I) 60 REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND) C 60 PRINT *,"DEBUG: gone 60, REI(",L,",",I,") =",REI(L,I) RETURN 65 DO 80 L=1,LI LPO=L+1 Loading Loading @@ -1138,7 +1093,6 @@ cccccccccccccccccccccc 80 DREF(L)=DY2*SZ+Y2 CRI=(1.0D0,0.0D0) IF(MOD(NSH,2).NE.0)CRI=DC0(IC)/EXDC C PRINT *,"DEBUG: going 90, AREX =",AREX DO 90 L=1,LI LPO=L+1 LTPO=LPO+L Loading @@ -1150,13 +1104,11 @@ C PRINT *,"DEBUG: going 90, AREX =",AREX CCNC=RMF(L)*DFB CCND=DRMF(L)*FB(LPO)*SZ*LTPO RMI(L,I)=1.0D0+UIM*(CCNA-CCNB)/(CCNC-CCND) C PRINT *,"DEBUG: gone 90, RMI(",L,",",I,") =",RMI(L,I) CCNA=REF(L)*DFN CCNB=DREF(L)*FN(LPO)*SZ*LTPO CCNC=REF(L)*DFB CCND=DREF(L)*FB(LPO)*SZ*LTPO 90 REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND) C 90 PRINT *,"DEBUG: gone 90, REI(",L,",",I,") =",REI(L,I) RETURN END SUBROUTINE RKT(NPNTMO,STEP,X,LPO,Y1,Y2,DY1,DY2) Loading Loading @@ -1498,18 +1450,13 @@ CCC LL=LM PI4=DACOS(0.0D0)*8.0D0 PI4IRS=1.0D0/DSQRT(PI4) X=COSRTH C PRINT *,"DEBUG: X =",X Y=DABS(SINRTH) C PRINT *,"DEBUG: Y =",Y CLLMO=3.0D0 CLL=1.5D0 YTOL=Y PLEGN(1)=1.0D0 C PRINT *,"DEBUG: PLEGN( 1 ) =",PLEGN(1) PLEGN(2)=X*DSQRT(CLLMO) C PRINT *,"DEBUG: PLEGN( 2 ) =",PLEGN(2) PLEGN(3)=YTOL*DSQRT(CLL) C PRINT *,"DEBUG: PLEGN( 3 ) =",PLEGN(3) SINRMP(1)=SINRPH COSRMP(1)=COSRPH IF(LL.LT.2)GO TO 30 Loading @@ -1529,17 +1476,14 @@ C PRINT *,"DEBUG: PLEGN( 3 ) =",PLEGN(3) CDM=(LTMO-2)*LS 10 PLEGN(MPOPK)=PLEGN(MPOPK-L)*X*DSQRT(CN/CD)- 1PLEGN(MPOPK-LTMO)*DSQRT(CNM/CDM) C10 PRINT *,"DEBUG: PLEGN(",MPOPK,") =",PLEGN(MPOPK) LPK=L+K CLTPO=LTPO PLEGN(LPK)=PLEGN(K)*X*DSQRT(CLTPO) C PRINT *,"DEBUG: PLEGN(",LPK,") =",PLEGN(LPK) K=LPK+1 CLT=LTPO-1 CLL=CLL*(CLTPO/CLT) YTOL=YTOL*Y PLEGN(K)=YTOL*DSQRT(CLL) C PRINT *,"DEBUG: PLEGN(",K,") =",PLEGN(K) SINRMP(L)=SINRPH*COSRMP(LMO)+COSRPH*SINRMP(LMO) 20 COSRMP(L)=COSRPH*COSRMP(LMO)-SINRPH*SINRMP(LMO) 30 L=0 Loading @@ -1548,17 +1492,14 @@ C PRINT *,"DEBUG: PLEGN(",K,") =",PLEGN(K) L0Y=K+1 L0P=K/2+1 YLM(L0Y)=PI4IRS*PLEGN(L0P) C PRINT *, "DEBUG: YLM(",L0Y,") =",YLM(L0Y) 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) C PRINT *, "DEBUG: YLM(",LMY,") =",YLM(LMY) LMY=L0Y-M YLM(LMY)=SAVE*DCMPLX(COSRMP(M),-SINRMP(M)) C PRINT *, "DEBUG: YLM(",LMY,") =",YLM(LMY) 45 IF(M.GE.L)GO TO 47 M=M+1 GO TO 44 Loading