diff --git a/cluster/clu.f b/cluster/clu.f
new file mode 100644
index 0000000000000000000000000000000000000000..8f1c2172f1534fff341322d0ab075e8689e61c74
--- /dev/null
+++ b/cluster/clu.f
@@ -0,0 +1,3469 @@
+ PROGRAM CLU
+CCC 111130
+CCC SUPPORTS EXPERIMENTAL DIELECTRIC FUNCTIONS ONLY
+CCC NSPH=6; LI=8, LE=8; NXI=200
+CCC
+CCC FOR IDFC>=0 ONLY,
+CCC WHEN JXI=JWTM WRITES THE JWTM-TH TRANSITION MATRIX
+CCC
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC COMMON/C1/VH(NCOU*LITPO),VJ0(NSPH*LMTPO),
+CCC 1VYHJ(NCOU*LITPOS),VYJ0(NSPH*LMTPOS),
+CCC 2RMI(LI,NSPH),REI(LI,NSPH),W(NLEMT,4),
+CCC 3AM0M(NLEMT,NLEMT),FSAS(NSPH),SAS(NSPH,2,2),VINTS(NSPH,16),
+CCC 4V3J0(NV3J),SSCS(NSPH),SEXS(NSPH),SABS(NSPH),
+CCC 5SQSCS(NSPH),SQEXS(NSPH),SQABS(NSPH),GCSV(NSPH),
+CCC 6RXX(NSPH),RYY(NSPH),RZZ(NSPH),ROS(NSPH),RC(NSPH,NSHL),
+CCC 7IND3J(LMPO,LM),IOG(NSPH),NSHL(NSPH)
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+CCC COMMON/C2/RIS(NHSPO),DLRI(NHSPO),DC0(NSHL-NTL+1),
+CCC 1VKT(NSPH),VSZ(NSPH)
+CCC NHSPO=MAX0(NPNT,NPNTTS)*2-1
+ COMMON/C2/RIS(999),DLRI(999),DC0(5),VKT(6),VSZ(6)
+ COMPLEX*16 RIS,DLRI,DC0,VKT
+ COMMON/C3/TFSAS,TSAS(2,2),GCS,SCS,ECS,ACS
+ COMPLEX*16 TFSAS,TSAS
+ COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+CCC COMMON/C6/RAC3J(LMTPO)
+ COMMON/C6/RAC3J(17)
+CCC COMMON/C9/GIS(NDI,NLEM),GLS(NDI,NLEM),SAM(NDIT,NLEMT)
+ COMMON/C9/GIS(480,80),GLS(480,80),SAM(960,160)
+ COMPLEX*16 GIS,GLS,SAM
+CCC DIMENSION AM(MXNDM,MXNDM)
+CCC NDIT.LE.MXNDM
+ DIMENSION AM(960,960)
+ COMPLEX*16 AM
+CCC DIMENSION RCF(NSPH,NSHL)
+ DIMENSION RCF(6,8),CMULLR(4,4),CMUL(4,4),CEXTLR(4,4),CEXT(4,4),
+ 1U(3),US(3),UN(3),UNS(3),UP(3),UPS(3),
+ 2UNMP(3),UNSMP(3),UPMP(3),UPSMP(3),DUK(3),ARGI(1),ARGS(1)
+CCC DIMENSION ZPV(LM,3,2,2),GAPS(NSPH)
+ DIMENSION ZPV(8,3,2,2),GAPS(6),GAPV(3),
+ 1GAP(3,2),GAPP(3,2),GAPM(3,2),GAPPM(3,2)
+ COMPLEX*16 GAPP,GAPPM
+CCC DIMENSION TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH)
+ DIMENSION TQSE(2,6),TQSPE(2,6),TQSS(2,6),TQSPS(2,6),
+ 1TQEV(3),TQSV(3),TQCE(2,3),TQCPE(2,3),TQCS(2,3),TQCPS(2,3)
+ COMPLEX*16 TQSPE,TQSPS,TQCPE,TQCPS
+CCC DIMENSION DC0M(NSPH,NSHL-NTL),XIV(NXI)
+ DIMENSION DC0M(6,4),XIV(200)
+ COMPLEX*16 DC0M
+ COMPLEX*16 ARG,S0,S0M
+ 8010 FORMAT(1H ,16I5)
+ 8015 FORMAT(1H ,
+ 1'READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM')
+ 8050 FORMAT(3(1PD17.8))
+ 8060 FORMAT(1H ,'READ(IR,*)RXX(I),RYY(I),RZZ(I)')
+ 8068 FORMAT(1H ,'READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST')
+ 8069 FORMAT(1H ,'READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST')
+ 8170 FORMAT(1H ,7(1PD10.3))
+ 8171 FORMAT(1H ,'READ(IR,*)JWTM')
+ 6000 FORMAT(1H )
+ 6050 FORMAT(1H ,' SIZE=',1PD15.7,', REFRACTIVE INDEX=',2(1PD15.7))
+ 6055 FORMAT(1H ,' SIZE=',1PD15.7)
+ 6060 FORMAT(1H ,' CIRC ',I2)
+ 6061 FORMAT(1H ,' CIRC')
+ 6065 FORMAT(1H ,' LIN ',I2)
+ 6066 FORMAT(1H ,' LIN')
+ 6068 FORMAT(1H ,4X,'SPHERES; LMX=LI')
+ 6069 FORMAT(1H ,4X,'SPHERES; LMX=MIN0(LI,LE)')
+ 6070 FORMAT(1H ,4X,'SPHERE ',I2)
+ 6080 FORMAT(1H ,4X,'CLUSTER')
+ 6100 FORMAT(1H ,'----- SCS ----- ABS ----- EXS ----- ALBEDS --')
+ 6150 FORMAT(1H ,1PD14.7,5(1PD15.7))
+ 6165 FORMAT(1H ,' UN=(',1PD12.5,',',1PD12.5,',',1PD12.5,')')
+ 6166 FORMAT(1H ,' UNI=(',1PD12.5,',',1PD12.5,',',1PD12.5,')')
+ 6167 FORMAT(1H ,' UNS=(',1PD12.5,',',1PD12.5,',',1PD12.5,')')
+ 6200 FORMAT(1H ,'---- SCS/GS -- ABS/GS -- EXS/GS ---')
+ 6240 FORMAT(1H ,1X,'FSAS=',2(1PD15.7))
+ 6242 FORMAT(1H ,1X,'SAS(1,1)=',2(1PD15.7),', SAS(2,1)=',2(1PD15.7))
+ 6243 FORMAT(1H ,1X,'SAS(1,2)=',2(1PD15.7),', SAS(2,2)=',2(1PD15.7))
+ 6250 FORMAT(1H ,1X,'FSAT=',2(1PD15.7))
+ 6252 FORMAT(1H ,1X,'SAT(1,1)=',2(1PD15.7),', SAT(2,1)=',2(1PD15.7))
+ 6253 FORMAT(1H ,1X,'SAT(1,2)=',2(1PD15.7),', SAT(2,2)=',2(1PD15.7))
+ 6300 FORMAT(1H ,' STOP IN LUCIN: MATRIX SINGULAR')
+ 6331 FORMAT(1H ,1X,'MULS')
+ 6332 FORMAT(1H ,1X,'MULC')
+ 6336 FORMAT(1H ,1X,'MULSLR')
+ 6337 FORMAT(1H ,1X,'MULCLR')
+ 6350 FORMAT((1H ,7X,4(1PD15.7)))
+ 6460 FORMAT(1H ,4X,'CLUSTER (ENSEMBLE AVERAGE, MODE',I2,')')
+ 6500 FORMAT(1H ,'----- SCC ----- ABC ----- EXC ----- ALBEDC --')
+ 6550 FORMAT
+ 1(1H ,' AT ',I1,' LCALC=',I3,' TOO SMALL WITH ARG=',2(1PD15.7))
+ 6600 FORMAT(1H ,'--- SCC/TGS - ABC/TGS - EXC/TGS ---')
+ 6650 FORMAT(1H ,' STOP IN HJV')
+ 6700 FORMAT(1H ,'---- SCCRT --- ABCRT --- EXCRT ----')
+ 6710 FORMAT(1H ,' REFR. INDEX OF EXTERNAL MEDIUM=',1PD15.7)
+ 6714 FORMAT(1H ,' VK=',1PD15.7,', XI IS SCALE FACTOR FOR LENGTHS')
+ 6715 FORMAT(1H ,' VK=',1PD15.7,', XI=',1PD15.7)
+ 6716 FORMAT(1H ,' XI=',1PD15.7)
+ 6720 FORMAT(1H ,' TIDG=',1PD10.3,', PIDG=',1PD10.3,
+ 1', TSDG=',1PD10.3,', PSDG=',1PD10.3)
+ 6721 FORMAT(1H ,' CFMP=',1PD15.7,', SFMP=',1PD15.7)
+ 6722 FORMAT(1H ,' CFSP=',1PD15.7,', SFSP=',1PD15.7)
+ 6740 FORMAT(1H ,' SCAND=',1PD10.3)
+ 6750 FORMAT(1H ,' STOP IN DME')
+ 6800 FORMAT(1H ,1X,2('FSAC(',I1,',',I1,')=',2(1PD15.7),3X))
+ 6810 FORMAT(1H ,1X,2(' SAC(',I1,',',I1,')=',2(1PD15.7),3X))
+ 6850 FORMAT(1H ,1X,'RE(FSAC(',I1,',',I1,'))/RE(TFSAS)=',1PD15.7,
+ 1', IM(FSAC(',I1,',',I1,'))/IM(TFSAS)=',1PD15.7)
+ 6930 FORMAT(1H ,1X,'QSCHU=',1PD15.7,', PSCHU=',1PD15.7,
+ 1', S0MAG=',1PD15.7)
+ 6940 FORMAT(1H ,1X,'(RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=',
+ 11PD15.7)
+ 6945 FORMAT(1H ,1X,'(IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=',
+ 11PD15.7)
+ 6970 FORMAT(1H ,' COSAV=',1PD15.7,', RAPRS=',1PD15.7)
+ 6971 FORMAT(1H ,' Fx=',1PD15.7,', Fy=',1PD15.7,', Fz=',1PD15.7)
+ 6973 FORMAT(1H ,' Fl=',1PD15.7,', Fr=',1PD15.7,', Fk=',1PD15.7)
+ 6974 FORMAT(1H ,' Fk=',1PD15.7)
+ 6977 FORMAT(1H ,' IPO=',I2,', TQEk=',1PD15.7,', TQSk=',1PD15.7)
+ 6981 FORMAT
+ 1(1H ,' TQEx=',1PD15.7,', TQEy=',1PD15.7,', TQEz=',1PD15.7)
+ 6982 FORMAT
+ 1(1H ,' TQSx=',1PD15.7,', TQSy=',1PD15.7,', TQSz=',1PD15.7)
+ 6983 FORMAT
+ 1(1H ,' TQEl=',1PD15.7,', TQEr=',1PD15.7,', TQEk=',1PD15.7)
+ 6984 FORMAT
+ 1(1H ,' TQSl=',1PD15.7,', TQSr=',1PD15.7,', TQSk=',1PD15.7)
+ 6991 FORMAT('========== JXI =',I3,' ====================')
+ 6992 FORMAT('********** JTH =',I3,', JPH =',I3,
+ 1', JTHS =',I3,', JPHS =',I3,' ********************')
+ IR=5
+ IW=6
+ IT=7
+ ITIN=17
+CCC
+CCC OTHER COMMENTS IN FILE GLOBALCOMS
+CCC
+CCC
+CCC CONSTANTS APPEAR AS (...).X(...)D0, OR (...).X(...)D+X(X),
+CCC OR (...).X(...)D-X(X)
+CCC
+CCC
+CCC MXNDM IS FIRST NUMERICAL DECLARED DIMENSION OF AM;
+CCC NDIT=2*NDI; NDI=NSPH*NLIM
+CCC
+CCC
+CCC VK*ROS(I) IS SIZE PARAMETER in vacuo
+CCC
+CCC
+CCC OUTPUTS (NORMALIZED SCATTERING AMPLITUDE)*(2*PIG/(VK*VK)) IN m**3;
+CCC OUTPUTS FORCE*(c/I0) IN m**2;
+CCC OUTPUTS (EXT_TO_TORQUE) AND (SCA_TO_TORQUE) CONTRIBUTIONS
+CCC *(c/I0)*(VK*EXRI)**3; (TORQUE OUTPUTS)/(EXDC*VK*VK) ARE IN m**2
+CCC
+CCC
+CCC ON A LINEAR BASIS,
+CCC POLARIZATION (1), (1,ni), (ns,1) OR -1 IS PARALLE(l),
+CCC POLARIZATION (2), (2,ni), (ns,2) OR 1 IS PERPENDICULA(r);
+CCC ON A CIRCULAR BASIS,
+CCC POLARIZATION (1), (1,ni), (ns,1) OR -1 IS (R)ight,
+CCC POLARIZATION (2), (2,ni), (ns,2) OR 1 IS (L)eft;
+CCC ns=SCATWAVE POLARIZATION, ni=INCWAVE POLARIZATION
+CCC
+CCC
+CCC NV3J=(LM*LMPO*(LMT+7))/6;
+CCC LM=MAX0(LI,LE)
+CCC
+CCC
+CCC IAVM=0 --> RANDOM AVERAGES FOR AMPLITUDE AND CROSS SECTIONS;
+CCC IAVM=1 --> SAME AS IAVM=0, PLUS RANDOM AVERAGES FOR INTENSITY
+CCC
+CCC
+CCC ISAM<0 --> ARBITRARILY VARYING SCATTERING PLANE,
+CCC AS IT IS UNDEFINED FOR FORWARD OR BACKWARD SCATTERING,
+CCC MERIDIONAL PLANE FOR INCIDENCE ANGLES IS ASSUMED;
+CCC ISAM=0 --> MERIDIONAL PLANES
+CCC (I. E. REFERENCE PLANES ACROSS Z AXIS)
+CCC AS DETERMINED BY INCIDENCE AND OBSERVATION ANGLES;
+CCC ISAM=1 --> SCATTERING PLANE ACROSS THE Z AXIS; SCATTERING
+CCC AS A FUNCTION OF INCIDENCE AND OF OBSERVATION ANGLES;
+CCC NEEDS ONLY TH, PH, THS; THSLST GT 180 IF NECESSARY;
+CCC ISAM>1 --> SCATTERING PLANE ACROSS THE Z AXIS; SCATTERING
+CCC AS A FUNCTION OF INCIDENCE ANGLES
+CCC FOR A FIXED SCATTERING ANGLE THSCA=THS-TH;
+CCC NEEDS ONLY TH, PH, THS
+CCC
+ OPEN(IR,FILE='DCLU',STATUS='OLD')
+ READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM
+ DO 101 I=1,NSPH
+ 101 READ(IR,*)RXX(I),RYY(I),RZZ(I)
+ READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST
+ READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST
+ READ(IR,*)JWTM
+ CLOSE(IR)
+ OPEN(IW,FILE='OCLU',STATUS='UNKNOWN')
+ WRITE(IW,8015)
+ WRITE(IW,8010)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM
+ WRITE(IW,8060)
+ WRITE(IW,8050)(RXX(I),RYY(I),RZZ(I),I=1,NSPH)
+ WRITE(IW,8068)
+ WRITE(IW,8170)TH,THSTP,THLST,THS,THSSTP,THSLST
+ WRITE(IW,8069)
+ WRITE(IW,8170)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST
+ WRITE(IW,8171)
+ WRITE(IW,8010)JWTM
+ OPEN(ITIN,FILE='TEDF',FORM='UNFORMATTED',STATUS='OLD')
+ READ(ITIN)NSPHT
+ IF(NSPHT.NE.NSPH)GO TO 497
+ READ(ITIN)(IOG(I),I=1,NSPH)
+ READ(ITIN)EXDC,WP,XIP,IDFC,NXI
+ READ(ITIN)(XIV(I),I=1,NXI)
+ DO 103 I=1,NSPH
+ IF(IOG(I).LT.I)GO TO 103
+ READ(ITIN)NSHL(I),ROS(I)
+ NSH=NSHL(I)
+ READ(ITIN)(RCF(I,NS),NS=1,NSH)
+ 103 CONTINUE
+ WRITE(IW,*)' READ(ITIN)NSPHT'
+ WRITE(IW,*)' READ(ITIN)(IOG(I),I=1,NSPH)'
+ WRITE(IW,*)' READ(ITIN)EXDC,WP,XIP,IDFC,NXI'
+ WRITE(IW,*)' READ(ITIN)(XIV(I),I=1,NXI)'
+ WRITE(IW,*)' READ(ITIN)NSHL(I),ROS(I)'
+ WRITE(IW,*)' READ(ITIN)(RCF(I,NS),NS=1,NSH)'
+ WRITE(IW,6000)
+ C0=0.0D0
+ SML=1.0D-3
+ IF(THSTP.EQ.C0)NTH=0
+ IF(THSTP.NE.C0)NTH=IDINT((THLST-TH)/THSTP+SML)
+ NTH=NTH+1
+ IF(PHSTP.EQ.C0)NPH=0
+ IF(PHSTP.NE.C0)NPH=IDINT((PHLST-PH)/PHSTP+SML)
+ NPH=NPH+1
+ IF(ISAM.LE.1)GO TO 110
+ NTHS=1
+ THSCA=THS-TH
+ GO TO 111
+ 110 IF(THSSTP.EQ.C0)NTHS=0
+ IF(THSSTP.NE.C0)NTHS=IDINT((THSLST-THS)/THSSTP+SML)
+ NTHS=NTHS+1
+ IF(ISAM.LT.1)GO TO 112
+ 111 NPHS=1
+ GO TO 113
+ 112 IF(PHSSTP.EQ.C0)NPHS=0
+ IF(PHSSTP.NE.C0)NPHS=IDINT((PHSLST-PHS)/PHSSTP+SML)
+ NPHS=NPHS+1
+ 113 NK=NTH*NPH
+ NKS=NTHS*NPHS
+ NKKS=NK*NKS
+ TH1=TH
+ PH1=PH
+ THS1=THS
+ PHS1=PHS
+ CALL STR(RCF)
+ LM=MAX0(LI,LE)
+ CALL THDPS(LM,ZPV)
+ EXRI=DSQRT(EXDC)
+ WRITE(IW,6710)EXRI
+ OPEN(IT,FILE='TPPOAN',FORM='UNFORMATTED',STATUS='UNKNOWN')
+ WRITE(IT)IAVM,ISAM,INPOL
+ WRITE(IT)NXI,NTH,NPH,NTHS,NPHS
+
+C
+ NLEMT=NLEM+NLEM
+C
+
+ WN=WP/3.0D08
+ SQSFI=1.0D0
+ IF(IDFC.GE.0)GO TO 118
+ VK=XIP*WN
+ WRITE(IW,6714)VK
+ WRITE(IW,6000)
+ 118 DO 488 JXI=1,NXI
+ JAW=1
+ WRITE(IW,6991)JXI
+ XI=XIV(JXI)
+ IF(IDFC.LT.0)GO TO 119
+ VK=XI*WN
+ VKARG=VK
+ WRITE(IW,6715)VK,XI
+ GO TO 120
+ 119 VKARG=XI*VK
+ SQSFI=1.0D0/(XI*XI)
+ WRITE(IW,6716)XI
+ 120 CALL HJV(EXRI,VKARG,JER,LCALC,ARG)
+ IF(JER.EQ.0)GO TO 122
+ WRITE(IW,6650)
+ GO TO 490
+ 122 DO 132 I=1,NSPH
+ IOGI=IOG(I)
+ IF(IOGI.EQ.I)GO TO 125
+ DO 123 L=1,LI
+ RMI(L,I)=RMI(L,IOGI)
+ 123 REI(L,I)=REI(L,IOGI)
+ GO TO 132
+ 125 NSH=NSHL(I)
+ ICI=(NSH+1)/2
+ IF(IDFC.NE.0)GO TO 127
+ READ(ITIN)(DC0(IC),IC=1,ICI)
+ GO TO 130
+ 127 IF(JXI.EQ.1)READ(ITIN)(DC0M(I,IC),IC=1,ICI)
+ DO 128 IC=1,ICI
+ 128 DC0(IC)=DC0M(I,IC)
+ 130 IF(MOD(NSH,2).EQ.0)DC0(ICI+1)=EXDC
+ CALL DME(LI,I,NPNT,NPNTTS,VKARG,EXDC,EXRI,JER,LCALC,ARG)
+ IF(JER.EQ.0)GO TO 132
+ WRITE(IW,6750)
+ GO TO 490
+ 132 CONTINUE
+ CALL CMS(AM)
+ NDIT=2*NSPH*NLIM
+ CALL LUCIN(AM,MXNDM,NDIT,JER)
+ IF(JER.EQ.1)GO TO 495
+ CALL ZTM(AM)
+
+C
+ IF(IDFC.LT.0)GO TO 156
+ IF(JXI.NE.JWTM)GO TO 156
+ IS=1
+ ITS=8
+ OPEN(ITS,FILE='TTMS',FORM='UNFORMATTED',STATUS='UNKNOWN')
+ WRITE(ITS)IS,LM
+ WRITE(ITS)VK,EXRI
+ WRITE(ITS)((AM0M(I,J),J=1,NLEMT),I=1,NLEMT)
+ CLOSE(ITS)
+ 156 CONTINUE
+C
+
+ IF(INPOL.NE.0)GO TO 158
+ WRITE(IW,6066)
+ GO TO 160
+ 158 WRITE(IW,6061)
+ 160 CS0=0.25D0*VK*VK*VK/DACOS(C0)
+ CALL SCR0(VK,EXRI)
+ 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)
+ IF(LI.NE.LE)WRITE(IW,6068)
+ DO 170 I=1,NSPH
+ IF(IOG(I).LT.I)GO TO 170
+ ALBEDS=SSCS(I)/SEXS(I)
+ SQSCS(I)=SQSCS(I)*SQSFI
+ SQABS(I)=SQABS(I)*SQSFI
+ SQEXS(I)=SQEXS(I)*SQSFI
+ WRITE(IW,6070)I
+ IF(NSHL(I).EQ.1)GO TO 162
+ WRITE(IW,6055)VSZ(I)
+ GO TO 164
+ 162 WRITE(IW,6050)VSZ(I),VKT(I)
+ 164 WRITE(IW,6100)
+ WRITE(IW,6150)SSCS(I),SABS(I),SEXS(I),ALBEDS
+ WRITE(IW,6200)
+ WRITE(IW,6150)SQSCS(I),SQABS(I),SQEXS(I)
+ WRITE(IW,6240)FSAS(I)
+ CSCH=2.0D0*VK*SQSFI/GCSV(I)
+ S0=FSAS(I)*EXRI
+ QSCHU=DIMAG(S0)*CSCH
+ PSCHU=DREAL(S0)*CSCH
+ S0MAG=CDABS(S0)*CS0
+ WRITE(IW,6930)QSCHU,PSCHU,S0MAG
+ RAPR=SEXS(I)-GAPS(I)
+ COSAV=GAPS(I)/SSCS(I)
+ WRITE(IW,6970)COSAV,RAPR
+ WRITE(IW,6977)1,TQSE(1,I),TQSS(1,I)
+ WRITE(IW,6977)2,TQSE(2,I),TQSS(2,I)
+ 170 CONTINUE
+ WRITE(IW,6250)TFSAS
+ CSCH=2.0D0*VK*SQSFI/GCS
+ S0=TFSAS*EXRI
+ QSCHU=DIMAG(S0)*CSCH
+ PSCHU=DREAL(S0)*CSCH
+ S0MAG=CDABS(S0)*CS0
+ WRITE(IW,6930)QSCHU,PSCHU,S0MAG
+ WRITE(IT)VK
+ CALL PCRSM0(VK,EXRI,INPOL)
+ CALL APCRA(ZPV,LE,AM0M,INPOL,SQK,GAPM,GAPPM)
+ TH=TH1
+ DO 486 JTH=1,NTH
+ PH=PH1
+ DO 484 JPH=1,NPH
+ IF((NK.EQ.1).AND.(JXI.GT.1))GO TO 180
+ CALL UPVMP(TH,PH,0,COST,SINT,COSP,SINP,U,UPMP,UNMP)
+ IF(ISAM.LT.0)GO TO 184
+ CALL WMAMP
+ 1(0,COST,SINT,COSP,SINP,INPOL,LE,0,NSPH,ARGI,U,UPMP,UNMP)
+ GO TO 182
+ 180 IF(ISAM.LT.0)GO TO 184
+ 182 CALL APC(ZPV,LE,AM0M,W,SQK,GAP,GAPP)
+ CALL RABA(LE,AM0M,W,TQCE,TQCPE,TQCS,TQCPS)
+ JW=1
+ 184 THSL=THS1
+ DO 482 JTHS=1,NTHS
+ THS=THSL
+ ICSPNV=0
+ IF(ISAM.GT.1)THS=TH+THSCA
+ IF(ISAM.LT.1)GO TO 186
+ PHSPH=C0
+ IF((THS.LT.C0).OR.(THS.GT.180.0D0))PHSPH=180.0D0
+ IF(THS.LT.C0)THS=-THS
+ IF(THS.GT.180.0D0)THS=360.0D0-THS
+ IF(PHSPH.NE.C0)ICSPNV=1
+ 186 PHS=PHS1
+ DO 480 JPHS=1,NPHS
+ IF(ISAM.LT.1)GO TO 188
+ PHS=PH+PHSPH
+ IF(PHS.GE.360.0D0)PHS=PHS-360.0D0
+ 188 IF((NKS.EQ.1).AND.((JXI.GT.1).OR.(JTH.GT.1).OR.(JPH.GT.1)))
+ 1GO TO 190
+ CALL UPVMP(THS,PHS,ICSPNV,COSTS,SINTS,COSPS,SINPS,US,UPSMP,UNSMP)
+ IF(ISAM.GE.0)CALL WMAMP
+ 1(2,COSTS,SINTS,COSPS,SINPS,INPOL,LE,0,NSPH,ARGS,US,UPSMP,UNSMP)
+ 190 IF((NKKS.EQ.1).AND.(JXI.GT.1))GO TO 194
+ CALL UPVSP(U,UPMP,UNMP,US,UPSMP,UNSMP,UP,UN,UPS,UNS,DUK,
+ 1ISQ,IBF,SCAN,CFMP,SFMP,CFSP,SFSP)
+ IF(ISAM.GE.0)GO TO 192
+ CALL WMASP(COST,SINT,COSP,SINP,COSTS,SINTS,COSPS,SINPS,
+ 1U,UP,UN,US,UPS,UNS,ISQ,IBF,INPOL,LE,0,NSPH,ARGI,ARGS)
+ GO TO 194
+ 192 DO 193 I=1,3
+ UP(I)=UPMP(I)
+ UN(I)=UNMP(I)
+ UPS(I)=UPSMP(I)
+ 193 UNS(I)=UNSMP(I)
+ 194 IF(IAVM.EQ.1)CALL CRSM1(VK,EXRI)
+ IF(ISAM.GE.0)GO TO 196
+ CALL APC(ZPV,LE,AM0M,W,SQK,GAP,GAPP)
+ CALL RABA(LE,AM0M,W,TQCE,TQCPE,TQCS,TQCPS)
+ JW=1
+ 196 WRITE(IT)TH,PH,THS,PHS
+ WRITE(IT)SCAN
+ IF(JAW.EQ.0)GO TO 212
+ JAW=0
+ CALL MEXTC(VK,EXRI,FSACM,CEXTLR,CEXT)
+ WRITE(IT)((CEXT(I,J),J=1,4),I=1,4)
+ WRITE(IT)(SCSCM(I),SCSCPM(I),ECSCM(I),ECSCPM(I),I=1,2)
+ WRITE(IT)((GAPM(I,J),GAPPM(I,J),J=1,2),I=1,3)
+ WRITE(IW,6460)IAVM
+ JLR=2
+ DO 210 ILR=1,2
+ IPOL=(-1)**ILR
+ IF(ILR.EQ.2)JLR=1
+ EXTSM=ECSCM(ILR)
+ QEXTM=EXTSM*SQSFI/GCS
+ EXTRM=EXTSM/ECS
+ SCASM=SCSCM(ILR)
+ ALBDM=SCASM/EXTSM
+ QSCAM=SCASM*SQSFI/GCS
+ SCARM=SCASM/SCS
+ ABSSM=EXTSM-SCASM
+ QABSM=ABSSM*SQSFI/GCS
+ IF(DABS(ACS/ECS).GT.1.0D-06)GO TO 202
+ ABSRM=1.0D0
+ GO TO 204
+ 202 ABSRM=ABSSM/ACS
+ 204 S0M=FSACM(ILR,ILR)*EXRI
+ QSCHUM=DIMAG(S0M)*CSCH
+ PSCHUM=DREAL(S0M)*CSCH
+ S0MAGM=CDABS(S0M)*CS0
+ RFINRM=DREAL(FSACM(ILR,ILR))/DREAL(TFSAS)
+ EXTCRM=DIMAG(FSACM(ILR,ILR))/DIMAG(TFSAS)
+ IF(INPOL.NE.0)GO TO 206
+ WRITE(IW,6065)IPOL
+ GO TO 208
+ 206 WRITE(IW,6060)IPOL
+ 208 WRITE(IW,6500)
+ WRITE(IW,6150)SCASM,ABSSM,EXTSM,ALBDM
+ WRITE(IW,6600)
+ WRITE(IW,6150)QSCAM,QABSM,QEXTM
+ WRITE(IW,6700)
+ WRITE(IW,6150)SCARM,ABSRM,EXTRM
+ WRITE(IW,6800)ILR,ILR,FSACM(ILR,ILR),JLR,ILR,FSACM(JLR,ILR)
+ WRITE(IW,6850)ILR,ILR,RFINRM,ILR,ILR,EXTCRM
+ WRITE(IW,6930)QSCHUM,PSCHUM,S0MAGM
+ RAPR=ECSCM(ILR)-GAPM(3,ILR)
+ COSAV=GAPM(3,ILR)/SCSCM(ILR)
+ FZ=RAPR
+ WRITE(IW,6970)COSAV,RAPR
+ WRITE(IW,6974)FZ
+ 210 CONTINUE
+ RMBRIF=(DREAL(FSACM(1,1))-DREAL(FSACM(2,2)))/DREAL(FSACM(1,1))
+ RMDCHR=(DIMAG(FSACM(1,1))-DIMAG(FSACM(2,2)))/DIMAG(FSACM(1,1))
+ WRITE(IW,6940)RMBRIF
+ WRITE(IW,6945)RMDCHR
+ 212 WRITE(IW,6992)JTH,JPH,JTHS,JPHS
+ WRITE(IW,6720)TH,PH,THS,PHS
+ WRITE(IW,6740)SCAN
+ WRITE(IW,6721)CFMP,SFMP
+ WRITE(IW,6722)CFSP,SFSP
+ IF(ISAM.LT.0)GO TO 214
+ WRITE(IW,6166)UN(1),UN(2),UN(3)
+ WRITE(IW,6167)UNS(1),UNS(2),UNS(3)
+ GO TO 220
+ 214 WRITE(IW,6165)UN(1),UN(2),UN(3)
+ 220 IF(INPOL.NE.0)GO TO 222
+ WRITE(IW,6066)
+ GO TO 224
+ 222 WRITE(IW,6061)
+ 224 CALL SCR2(VK,VKARG,EXRI,DUK)
+ IF(LI.NE.LE)WRITE(IW,6069)
+ DO 226 I=1,NSPH
+ IF(IOG(I).LT.I)GO TO 226
+ WRITE(IW,6070)I
+ WRITE(IW,6242)SAS(I,1,1),SAS(I,2,1)
+ WRITE(IW,6243)SAS(I,1,2),SAS(I,2,2)
+ DO 225 J=1,16
+ 225 VINT(J)=VINTS(I,J)
+ CALL MMULC(VINT,CMULLR,CMUL)
+ WRITE(IW,6331)
+ WRITE(IW,6350)((CMUL(I1,J),J=1,4),I1=1,4)
+ WRITE(IW,6336)
+ WRITE(IW,6350)((CMULLR(I1,J),J=1,4),I1=1,4)
+ 226 CONTINUE
+ WRITE(IW,6252)TSAS(1,1),TSAS(2,1)
+ WRITE(IW,6253)TSAS(1,2),TSAS(2,2)
+ WRITE(IW,6080)
+ CALL PCROS(VK,EXRI)
+ CALL MEXTC(VK,EXRI,FSAC,CEXTLR,CEXT)
+ CALL MMULC(VINT,CMULLR,CMUL)
+ IF(JW.EQ.0)GO TO 254
+ JW=0
+ WRITE(IT)((CEXT(I,J),J=1,4),I=1,4)
+ WRITE(IT)(SCSC(I),SCSCP(I),ECSC(I),ECSCP(I),I=1,2)
+ WRITE(IT)((GAP(I,J),GAPP(I,J),J=1,2),I=1,3)
+ WRITE(IT)((TQCE(I,J),TQCPE(I,J),J=1,3),I=1,2)
+ WRITE(IT)((TQCS(I,J),TQCPS(I,J),J=1,3),I=1,2)
+ WRITE(IT)(U(I),UP(I),UN(I),I=1,3)
+ 254 WRITE(IT)(VINT(I),I=1,16)
+ WRITE(IT)((CMUL(I,J),J=1,4),I=1,4)
+ JLR=2
+ DO 290 ILR=1,2
+ IPOL=(-1)**ILR
+ IF(ILR.EQ.2)JLR=1
+ EXTSEC=ECSC(ILR)
+ QEXT=EXTSEC*SQSFI/GCS
+ EXTRAT=EXTSEC/ECS
+ SCASEC=SCSC(ILR)
+ ALBEDC=SCASEC/EXTSEC
+ QSCA=SCASEC*SQSFI/GCS
+ SCARAT=SCASEC/SCS
+ ABSSEC=EXTSEC-SCASEC
+ QABS=ABSSEC*SQSFI/GCS
+ IF(DABS(ACS/ECS).GT.1.0D-06)GO TO 260
+ ABSRAT=1.0D0
+ GO TO 262
+ 260 ABSRAT=ABSSEC/ACS
+ 262 S0=FSAC(ILR,ILR)*EXRI
+ QSCHU=DIMAG(S0)*CSCH
+ PSCHU=DREAL(S0)*CSCH
+ S0MAG=CDABS(S0)*CS0
+ REFINR=DREAL(FSAC(ILR,ILR))/DREAL(TFSAS)
+ EXTCOR=DIMAG(FSAC(ILR,ILR))/DIMAG(TFSAS)
+ IF(INPOL.NE.0)GO TO 273
+ WRITE(IW,6065)IPOL
+ GO TO 275
+ 273 WRITE(IW,6060)IPOL
+ 275 WRITE(IW,6500)
+ WRITE(IW,6150)SCASEC,ABSSEC,EXTSEC,ALBEDC
+ WRITE(IW,6600)
+ WRITE(IW,6150)QSCA,QABS,QEXT
+ WRITE(IW,6700)
+ WRITE(IW,6150)SCARAT,ABSRAT,EXTRAT
+ WRITE(IW,6800)ILR,ILR,FSAC(ILR,ILR),JLR,ILR,FSAC(JLR,ILR)
+ WRITE(IW,6810)ILR,ILR,SAC(ILR,ILR),JLR,ILR,SAC(JLR,ILR)
+ WRITE(IW,6850)ILR,ILR,REFINR,ILR,ILR,EXTCOR
+ WRITE(IW,6930)QSCHU,PSCHU,S0MAG
+ IF((ISAM.GE.0).AND.((JTHS.GT.1).OR.(JPHS.GT.1)))GO TO 290
+ GAPV(1)=GAP(1,ILR)
+ GAPV(2)=GAP(2,ILR)
+ GAPV(3)=GAP(3,ILR)
+ EXTINS=ECSC(ILR)
+ SCATTS=SCSC(ILR)
+ CALL RFTR
+ 1(U,UP,UN,GAPV,EXTINS,SCATTS,RAPR,COSAV,FP,FN,FK,FX,FY,FZ)
+ WRITE(IW,6970)COSAV,RAPR
+ WRITE(IW,6973)FP,FN,FK
+ WRITE(IW,6971)FX,FY,FZ
+ TQEV(1)=TQCE(ILR,1)
+ TQEV(2)=TQCE(ILR,2)
+ TQEV(3)=TQCE(ILR,3)
+ TQSV(1)=TQCS(ILR,1)
+ TQSV(2)=TQCS(ILR,2)
+ TQSV(3)=TQCS(ILR,3)
+ CALL TQR(U,UP,UN,TQEV,TQSV,TEP,TEN,TEK,TSP,TSN,TSK)
+ WRITE(IW,6983)TEP,TEN,TEK
+ WRITE(IW,6984)TSP,TSN,TSK
+ WRITE(IW,6981)TQCE(ILR,1),TQCE(ILR,2),TQCE(ILR,3)
+ WRITE(IW,6982)TQCS(ILR,1),TQCS(ILR,2),TQCS(ILR,3)
+ 290 CONTINUE
+ RBIRIF=(DREAL(FSAC(1,1))-DREAL(FSAC(2,2)))/DREAL(FSAC(1,1))
+ RDICHR=(DIMAG(FSAC(1,1))-DIMAG(FSAC(2,2)))/DIMAG(FSAC(1,1))
+ WRITE(IW,6940)RBIRIF
+ WRITE(IW,6945)RDICHR
+ WRITE(IW,6332)
+ WRITE(IW,6350)((CMUL(I,J),J=1,4),I=1,4)
+ WRITE(IW,6337)
+ WRITE(IW,6350)((CMULLR(I,J),J=1,4),I=1,4)
+ IF(IAVM.EQ.0)GO TO 420
+ CALL MMULC(VINTM,CMULLR,CMUL)
+ WRITE(IT)(VINTM(I),I=1,16)
+ WRITE(IT)((CMUL(I,J),J=1,4),I=1,4)
+ WRITE(IW,6460)IAVM
+ IF(INPOL.NE.0)GO TO 316
+ WRITE(IW,6066)
+ GO TO 318
+ 316 WRITE(IW,6061)
+ 318 WRITE(IW,6332)
+ WRITE(IW,6350)((CMUL(I,J),J=1,4),I=1,4)
+ WRITE(IW,6337)
+ WRITE(IW,6350)((CMULLR(I,J),J=1,4),I=1,4)
+ 420 CONTINUE
+ 480 IF(ISAM.LT.1)PHS=PHS+PHSSTP
+ 482 IF(ISAM.LE.1)THSL=THSL+THSSTP
+ 484 PH=PH+PHSTP
+ 486 TH=TH+THSTP
+ 488 CONTINUE
+ GO TO 500
+ 490 WRITE(IW,6550)JER,LCALC,ARG
+ GO TO 500
+ 495 WRITE(IW,6300)
+ GO TO 500
+ 497 WRITE(IW,*)' WRONG INPUT TAPE'
+ 500 CONTINUE
+ CLOSE(IW)
+ CLOSE(IT)
+ CLOSE(ITIN)
+ STOP
+ END
+ SUBROUTINE RABAS(INPOL,LI,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS)
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC DIMENSION IOG(NSPH),RMI(LI,NSPH),REI(LI,NSPH),
+CCC 1TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH)
+ DIMENSION IOG(6),RMI(8,6),REI(8,6),
+ 1TQSE(2,6),TQSPE(2,6),TQSS(2,6),TQSPS(2,6)
+ COMPLEX*16 RMI,REI,TQSPE,TQSPS
+ COMPLEX*16 CC0,UIM,RM,RE,PCE,PCS
+ C0=0.0D0
+ PIG2=DACOS(C0)*4.0D0
+ CC0=(0.0D0,0.0D0)
+ UIM=(0.0D0,1.0D0)
+ DO 80 I=1,NSPH
+ IF(IOG(I).LT.I)GO TO 80
+ TQSE(1,I)=C0
+ TQSE(2,I)=C0
+ TQSPE(1,I)=CC0
+ TQSPE(2,I)=CC0
+ TQSS(1,I)=C0
+ TQSS(2,I)=C0
+ TQSPS(1,I)=CC0
+ TQSPS(2,I)=CC0
+ DO 70 L=1,LI
+ FL=L+L+1
+ RM=1.0D0/RMI(L,I)
+ RMM=DREAL(RM*DCONJG(RM))
+ RE=1.0D0/REI(L,I)
+ REM=DREAL(RE*DCONJG(RE))
+ IF(INPOL.NE.0)GO TO 40
+ PCE=((RM+RE)*UIM)*FL
+ PCS=((RMM+REM)*FL)*UIM
+ TQSPE(1,I)=TQSPE(1,I)-PCE
+ TQSPE(2,I)=TQSPE(2,I)+PCE
+ TQSPS(1,I)=TQSPS(1,I)-PCS
+ TQSPS(2,I)=TQSPS(2,I)+PCS
+ GO TO 70
+ 40 CE=DREAL(RM+RE)*FL
+ CS=(RMM+REM)*FL
+ TQSE(1,I)=TQSE(1,I)-CE
+ TQSE(2,I)=TQSE(2,I)+CE
+ TQSS(1,I)=TQSS(1,I)-CS
+ TQSS(2,I)=TQSS(2,I)+CS
+ 70 CONTINUE
+ IF(INPOL.NE.0)GO TO 75
+ TQSPE(1,I)=TQSPE(1,I)*PIG2
+ TQSPE(2,I)=TQSPE(2,I)*PIG2
+ TQSPS(1,I)=TQSPS(1,I)*PIG2
+ TQSPS(2,I)=TQSPS(2,I)*PIG2
+ 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
+ 80 CONTINUE
+ RETURN
+ END
+ SUBROUTINE TQR(U,UP,UN,TQEV,TQSV,TEP,TEN,TEK,TSP,TSN,TSK)
+ IMPLICIT REAL*8(A-H,O-Z)
+ DIMENSION U(3),UP(3),UN(3),TQEV(3),TQSV(3)
+ TEP=UP(1)*TQEV(1)+UP(2)*TQEV(2)+UP(3)*TQEV(3)
+ TEN=UN(1)*TQEV(1)+UN(2)*TQEV(2)+UN(3)*TQEV(3)
+ TEK=U(1)*TQEV(1)+U(2)*TQEV(2)+U(3)*TQEV(3)
+ TSP=UP(1)*TQSV(1)+UP(2)*TQSV(2)+UP(3)*TQSV(3)
+ TSN=UN(1)*TQSV(1)+UN(2)*TQSV(2)+UN(3)*TQSV(3)
+ TSK=U(1)*TQSV(1)+U(2)*TQSV(2)+U(3)*TQSV(3)
+ RETURN
+ END
+ SUBROUTINE RABA(LE,AM0M,W,TQCE,TQCPE,TQCS,TQCPS)
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC DIMENSION AM0M(NLEMT,NLEMT),W(NLEMT,4)
+ DIMENSION AM0M(160,160),W(160,4),
+ 1TQCE(2,3),TQCPE(2,3),TQCS(2,3),TQCPS(2,3)
+ COMPLEX*16 AM0M,W,TQCPE,TQCPS
+CCC DIMENSION A(NLEMT,2)
+ DIMENSION A(160,2),CTQCE(2,3),CTQCS(2,3)
+ COMPLEX*16 A,CTQCE,CTQCS
+ COMPLEX*16 CC0,UIM,ACW,ACWP,ACA,ACAP,C1,C2,C3
+ CC0=(0.0D0,0.0D0)
+ UIM=(0.0D0,1.0D0)
+ SQ2I=1.0D0/DSQRT(2.0D0)
+ NLEM=LE*(LE+2)
+ NLEMT=NLEM+NLEM
+ DO 20 I=1,NLEMT
+ C1=CC0
+ C2=CC0
+ DO 10 J=1,NLEMT
+ C1=C1+AM0M(I,J)*W(J,1)
+ C2=C2+AM0M(I,J)*W(J,2)
+ 10 CONTINUE
+ A(I,1)=C1
+ A(I,2)=C2
+ 20 CONTINUE
+ JPO=2
+ DO 70 IPO=1,2
+ IF(IPO.EQ.2)JPO=1
+ CTQCE(IPO,1)=CC0
+ CTQCE(IPO,2)=CC0
+ CTQCE(IPO,3)=CC0
+ TQCPE(IPO,1)=CC0
+ TQCPE(IPO,2)=CC0
+ TQCPE(IPO,3)=CC0
+ CTQCS(IPO,1)=CC0
+ CTQCS(IPO,2)=CC0
+ CTQCS(IPO,3)=CC0
+ TQCPS(IPO,1)=CC0
+ TQCPS(IPO,2)=CC0
+ TQCPS(IPO,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(A(I,IPO))*W(IMMU,IPO)+DCONJG(A(IE,IPO))*W(IMMUE,IPO)
+ ACWP=DCONJG(A(I,IPO))*W(IMMU,JPO)+DCONJG(A(IE,IPO))*W(IMMUE,JPO)
+ ACA=DCONJG(A(I,IPO))*A(IMMU,IPO)+DCONJG(A(IE,IPO))*A(IMMUE,IPO)
+ ACAP=DCONJG(A(I,IPO))*A(IMMU,JPO)+DCONJG(A(IE,IPO))*A(IMMUE,JPO)
+ CTQCE(IPO,1)=CTQCE(IPO,1)+ACW*RMU
+ TQCPE(IPO,1)=TQCPE(IPO,1)+ACWP*RMU
+ CTQCS(IPO,1)=CTQCS(IPO,1)+ACA*RMU
+ TQCPS(IPO,1)=TQCPS(IPO,1)+ACAP*RMU
+ 30 RMU=-M
+ ACW=DCONJG(A(I,IPO))*W(I,IPO)+DCONJG(A(IE,IPO))*W(IE,IPO)
+ ACWP=DCONJG(A(I,IPO))*W(I,JPO)+DCONJG(A(IE,IPO))*W(IE,JPO)
+ ACA=DCONJG(A(I,IPO))*A(I,IPO)+DCONJG(A(IE,IPO))*A(IE,IPO)
+ ACAP=DCONJG(A(I,IPO))*A(I,JPO)+DCONJG(A(IE,IPO))*A(IE,JPO)
+ CTQCE(IPO,2)=CTQCE(IPO,2)+ACW*RMU
+ TQCPE(IPO,2)=TQCPE(IPO,2)+ACWP*RMU
+ CTQCS(IPO,2)=CTQCS(IPO,2)+ACA*RMU
+ TQCPS(IPO,2)=TQCPS(IPO,2)+ACAP*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(A(I,IPO))*W(IMMU,IPO)+DCONJG(A(IE,IPO))*W(IMMUE,IPO)
+ ACWP=DCONJG(A(I,IPO))*W(IMMU,JPO)+DCONJG(A(IE,IPO))*W(IMMUE,JPO)
+ ACA=DCONJG(A(I,IPO))*A(IMMU,IPO)+DCONJG(A(IE,IPO))*A(IMMUE,IPO)
+ ACAP=DCONJG(A(I,IPO))*A(IMMU,JPO)+DCONJG(A(IE,IPO))*A(IMMUE,JPO)
+ CTQCE(IPO,3)=CTQCE(IPO,3)+ACW*RMU
+ TQCPE(IPO,3)=TQCPE(IPO,3)+ACWP*RMU
+ CTQCS(IPO,3)=CTQCS(IPO,3)+ACA*RMU
+ TQCPS(IPO,3)=TQCPS(IPO,3)+ACAP*RMU
+ 50 CONTINUE
+ 60 CONTINUE
+ 70 CONTINUE
+ DO 78 IPO=1,2
+ TQCE(IPO,1)=DREAL(CTQCE(IPO,1)-CTQCE(IPO,3))*SQ2I
+ TQCE(IPO,2)=DREAL((CTQCE(IPO,1)+CTQCE(IPO,3))*UIM)*SQ2I
+ TQCE(IPO,3)=DREAL(CTQCE(IPO,2))
+ C1=TQCPE(IPO,1)
+ C2=TQCPE(IPO,2)
+ C3=TQCPE(IPO,3)
+ TQCPE(IPO,1)=(C1-C3)*SQ2I
+ TQCPE(IPO,2)=(C1+C3)*(UIM*SQ2I)
+ TQCPE(IPO,3)=C2
+ TQCS(IPO,1)=-DREAL(CTQCS(IPO,1)-CTQCS(IPO,3))*SQ2I
+ TQCS(IPO,2)=-DREAL((CTQCS(IPO,1)+CTQCS(IPO,3))*UIM)*SQ2I
+ TQCS(IPO,3)=-DREAL(CTQCS(IPO,2))
+ C1=TQCPS(IPO,1)
+ C2=TQCPS(IPO,2)
+ C3=TQCPS(IPO,3)
+ TQCPS(IPO,1)=-(C1-C3)*SQ2I
+ TQCPS(IPO,2)=-(C1+C3)*(UIM*SQ2I)
+ TQCPS(IPO,3)=-C2
+ 78 CONTINUE
+ RETURN
+ END
+ SUBROUTINE RFTR
+ 1(U,UP,UN,GAPV,EXTINS,SCATTS,RAPR,COSAV,FP,FN,FK,FX,FY,FZ)
+ IMPLICIT REAL*8(A-H,O-Z)
+ DIMENSION U(3),UP(3),UN(3),GAPV(3)
+ FK=U(1)*GAPV(1)+U(2)*GAPV(2)+U(3)*GAPV(3)
+ RAPR=EXTINS-FK
+ COSAV=FK/SCATTS
+ FP=-(UP(1)*GAPV(1)+UP(2)*GAPV(2)+UP(3)*GAPV(3))
+ FN=-(UN(1)*GAPV(1)+UN(2)*GAPV(2)+UN(3)*GAPV(3))
+ FK=RAPR
+ FX=U(1)*EXTINS-GAPV(1)
+ FY=U(2)*EXTINS-GAPV(2)
+ FZ=U(3)*EXTINS-GAPV(3)
+ RETURN
+ END
+ SUBROUTINE MEXTC(VK,EXRI,FSAC,CEXTLR,CEXT)
+ IMPLICIT REAL*8(A-H,O-Z)
+ DIMENSION FSAC(2,2),CEXTLR(4,4),CEXT(4,4)
+ COMPLEX*16 FSAC
+ FA11R=DREAL(FSAC(1,1))
+ FA11I=DIMAG(FSAC(1,1))
+ FA21R=DREAL(FSAC(2,1))
+ FA21I=DIMAG(FSAC(2,1))
+ FA12R=DREAL(FSAC(1,2))
+ FA12I=DIMAG(FSAC(1,2))
+ FA22R=DREAL(FSAC(2,2))
+ FA22I=DIMAG(FSAC(2,2))
+ CEXTLR(1,1)=FA11I*2.0D0
+ CEXTLR(1,2)=0.0D0
+ CEXTLR(1,3)=-FA12I
+ CEXTLR(1,4)=-FA12R
+ CEXTLR(2,1)=0.0D0
+ CEXTLR(2,2)=FA22I*2.0D0
+ CEXTLR(2,3)=-FA21I
+ CEXTLR(2,4)=FA21R
+ CEXTLR(3,1)=-FA21I*2.0D0
+ CEXTLR(3,2)=-FA12I*2.0D0
+ CEXTLR(3,3)=FA11I+FA22I
+ CEXTLR(3,4)=FA22R-FA11R
+ CEXTLR(4,1)=FA21R*2.0D0
+ CEXTLR(4,2)=-FA12R*2.0D0
+ CEXTLR(4,3)=FA11R-FA22R
+ CEXTLR(4,4)=CEXTLR(3,3)
+ CEXT(1,1)=CEXTLR(4,4)
+ CEXT(2,2)=CEXTLR(4,4)
+ CEXT(3,3)=CEXTLR(4,4)
+ CEXT(3,4)=CEXTLR(3,4)
+ CEXT(4,3)=CEXTLR(4,3)
+ CEXT(4,4)=CEXTLR(4,4)
+ CEXT(1,2)=FA11I-FA22I
+ CEXT(1,3)=-FA12I-FA21I
+ CEXT(1,4)=FA21R-FA12R
+ CEXT(2,1)=CEXT(1,2)
+ CEXT(2,3)=FA21I-FA12I
+ CEXT(4,2)=FA12R+FA21R
+ CEXT(2,4)=-CEXT(4,2)
+ CEXT(3,1)=CEXT(1,3)
+ CEXT(3,2)=-CEXT(2,3)
+ CEXT(4,1)=CEXT(1,4)
+ CKM=VK/EXRI
+ DO 10 I=1,4
+ DO 10 J=1,4
+ CEXTLR(I,J)=CEXTLR(I,J)*CKM
+ 10 CEXT(I,J)=CEXT(I,J)*CKM
+ RETURN
+ END
+ SUBROUTINE MMULC(VINT,CMULLR,CMUL)
+ IMPLICIT REAL*8(A-H,O-Z)
+ DIMENSION VINT(16),CMULLR(4,4),CMUL(4,4)
+ COMPLEX*16 VINT
+ SM2=DREAL(VINT(1))
+ S24=DREAL(VINT(2))
+ D24=DIMAG(VINT(2))
+ SM4=DREAL(VINT(6))
+ S23=DREAL(VINT(9))
+ D32=DIMAG(VINT(9))
+ S34=DREAL(VINT(10))
+ D34=DIMAG(VINT(10))
+ SM3=DREAL(VINT(11))
+ S31=DREAL(VINT(12))
+ D31=DIMAG(VINT(12))
+ S21=DREAL(VINT(13))
+ D12=DIMAG(VINT(13))
+ S41=DREAL(VINT(14))
+ D14=DIMAG(VINT(14))
+ SM1=DREAL(VINT(16))
+ CMULLR(1,1)=SM2
+ CMULLR(1,2)=SM3
+ CMULLR(1,3)=-S23
+ CMULLR(1,4)=-D32
+ CMULLR(2,1)=SM4
+ CMULLR(2,2)=SM1
+ CMULLR(2,3)=-S41
+ CMULLR(2,4)=-D14
+ CMULLR(3,1)=-S24*2.0D0
+ CMULLR(3,2)=-S31*2.0D0
+ CMULLR(3,3)=S21+S34
+ CMULLR(3,4)=D34+D12
+ CMULLR(4,1)=-D24*2.0D0
+ CMULLR(4,2)=-D31*2.0D0
+ CMULLR(4,3)=D34-D12
+ CMULLR(4,4)=S21-S34
+ CMUL(1,1)=(SM2+SM3+SM4+SM1)*0.5D0
+ CMUL(1,2)=(SM2-SM3+SM4-SM1)*0.5D0
+ CMUL(1,3)=-S23-S41
+ CMUL(1,4)=-D32-D14
+ CMUL(2,1)=(SM2+SM3-SM4-SM1)*0.5D0
+ CMUL(2,2)=(SM2-SM3-SM4+SM1)*0.5D0
+ CMUL(2,3)=-S23+S41
+ CMUL(2,4)=-D32+D14
+ CMUL(3,1)=-S24-S31
+ CMUL(3,2)=-S24+S31
+ CMUL(3,3)=S21+S34
+ CMUL(3,4)=D34+D12
+ CMUL(4,1)=-D24-D31
+ CMUL(4,2)=-D24+D31
+ CMUL(4,3)=D34-D12
+ CMUL(4,4)=S21-S34
+ RETURN
+ END
+ SUBROUTINE SCR0(VK,EXRI)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ COMMON/C3/TFSAS,TSAS(2,2),GCS,SCS,ECS,ACS
+ COMPLEX*16 TFSAS,TSAS
+ COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+ COMPLEX*16 SUM21,RM,RE,CC0,CSAM
+ CC0=(0.0D0,0.0D0)
+ C0=0.0D0
+ EXDC=EXRI*EXRI
+ CCS=4.0D0*DACOS(C0)/(VK*VK)
+ CCCS=CCS/EXDC
+ CSAM=-(CCS/(EXRI*VK))*(0.0D0,0.5D0)
+ SCS=C0
+ ECS=C0
+ ACS=C0
+ TFSAS=CC0
+ DO 14 I=1,NSPH
+ IOGI=IOG(I)
+ IF(IOGI.LT.I)GO TO 12
+ SUMS=C0
+ SUM21=CC0
+ DO 10 L=1,LI
+ FL=L+L+1
+ RM=1.0D0/RMI(L,I)
+ RE=1.0D0/REI(L,I)
+ SUMS=SUMS+DREAL(DCONJG(RM)*RM+DCONJG(RE)*RE)*FL
+ 10 SUM21=SUM21+(RM+RE)*FL
+ SUM21=-SUM21
+ SCASEC=CCCS*SUMS
+ EXTSEC=-CCCS*DREAL(SUM21)
+ ABSSEC=EXTSEC-SCASEC
+ SSCS(I)=SCASEC
+ SEXS(I)=EXTSEC
+ SABS(I)=ABSSEC
+ GCSS=GCSV(I)
+ SQSCS(I)=SCASEC/GCSS
+ SQEXS(I)=EXTSEC/GCSS
+ SQABS(I)=ABSSEC/GCSS
+ FSAS(I)=SUM21*CSAM
+ 12 SCS=SCS+SSCS(IOGI)
+ ECS=ECS+SEXS(IOGI)
+ ACS=ACS+SABS(IOGI)
+ 14 TFSAS=TFSAS+FSAS(IOGI)
+ RETURN
+ END
+ SUBROUTINE SCR2(VK,VKARG,EXRI,DUK)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ COMMON/C3/TFSAS,TSAS(2,2),GCS,SCS,ECS,ACS
+ COMPLEX*16 TFSAS,TSAS
+ COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+ DIMENSION DUK(3)
+ COMPLEX*16 S11,S21,S12,S22,RM,RE,CC0,CSAM,CPH,PHAS,CC
+ CC0=(0.0D0,0.0D0)
+ C0=0.0D0
+ CCS=1.0D0/(VK*VK)
+ CSAM=-(CCS/(EXRI*VK))*(0.0D0,0.5D0)
+ PIGFSQ=6.4D+1*DACOS(C0)**2
+ CFSQ=4.0D0/(PIGFSQ*CCS*CCS)
+ CPH=(0.0D0,1.0D0)*(EXRI*VKARG)
+ LS=MIN0(LI,LE)
+ TSAS(1,1)=CC0
+ TSAS(2,1)=CC0
+ TSAS(1,2)=CC0
+ TSAS(2,2)=CC0
+ DO 14 I=1,NSPH
+ IOGI=IOG(I)
+ IF(IOGI.LT.I)GO TO 12
+ K=0
+ S11=CC0
+ S21=CC0
+ S12=CC0
+ S22=CC0
+ DO 10 L=1,LS
+ RM=1.0D0/RMI(L,I)
+ RE=1.0D0/REI(L,I)
+ LTPO=L+L+1
+ DO 10 IM=1,LTPO
+ K=K+1
+ KE=K+NLEM
+ S11=S11-W(K,3)*W(K,1)*RM-W(KE,3)*W(KE,1)*RE
+ S21=S21-W(K,4)*W(K,1)*RM-W(KE,4)*W(KE,1)*RE
+ 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
+ SAS(I,1,1)=S11*CSAM
+ SAS(I,2,1)=S21*CSAM
+ SAS(I,1,2)=S12*CSAM
+ SAS(I,2,2)=S22*CSAM
+ 12 PHAS=CDEXP(CPH*(DUK(1)*RXX(I)+DUK(2)*RYY(I)+DUK(3)*RZZ(I)))
+ TSAS(1,1)=TSAS(1,1)+SAS(IOGI,1,1)*PHAS
+ TSAS(2,1)=TSAS(2,1)+SAS(IOGI,2,1)*PHAS
+ TSAS(1,2)=TSAS(1,2)+SAS(IOGI,1,2)*PHAS
+ 14 TSAS(2,2)=TSAS(2,2)+SAS(IOGI,2,2)*PHAS
+ DO 24 I=1,NSPH
+ IOGI=IOG(I)
+ IF(IOGI.LT.I)GO TO 24
+ J=0
+ DO 22 IPO1=1,2
+ DO 22 JPO1=1,2
+ CC=DCONJG(SAS(I,JPO1,IPO1))
+ DO 22 IPO2=1,2
+ DO 22 JPO2=1,2
+ J=J+1
+ 22 VINTS(I,J)=SAS(I,JPO2,IPO2)*CC*CFSQ
+ 24 CONTINUE
+ J=0
+ DO 32 IPO1=1,2
+ DO 32 JPO1=1,2
+ CC=DCONJG(TSAS(JPO1,IPO1))
+ DO 32 IPO2=1,2
+ DO 32 JPO2=1,2
+ J=J+1
+ 32 VINTT(J)=TSAS(JPO2,IPO2)*CC*CFSQ
+ RETURN
+ END
+ SUBROUTINE PCROS(VK,EXRI)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+ COMPLEX*16 SUMP,SUM1,SUM2,SUM3,SUM4,AM,AMP,CC,CSAM,CC0
+ CC0=(0.0D0,0.0D0)
+ C0=0.0D0
+ EXDC=EXRI*EXRI
+ CCS=1.0D0/(VK*VK)
+ CCCS=CCS/EXDC
+ CSAM=-(CCS/(EXRI*VK))*(0.0D0,0.5D0)
+ PIGFSQ=6.4D+1*DACOS(C0)**2
+ CFSQ=4.0D0/(PIGFSQ*CCS*CCS)
+ NLEMT=NLEM+NLEM
+ JPO=2
+ DO 18 IPO=1,2
+ IF(IPO.EQ.2)JPO=1
+ IPOPT=IPO+2
+ JPOPT=JPO+2
+ SUM=C0
+ SUMP=CC0
+ SUM1=CC0
+ SUM2=CC0
+ SUM3=CC0
+ SUM4=CC0
+ DO 12 I=1,NLEMT
+ AM=CC0
+ AMP=CC0
+ DO 10 J=1,NLEMT
+ AM=AM+AM0M(I,J)*W(J,IPO)
+ 10 AMP=AMP+AM0M(I,J)*W(J,JPO)
+ SUM=SUM+DREAL(DCONJG(AM)*AM)
+ SUMP=SUMP+(DCONJG(AMP)*AM)
+ SUM1=SUM1+DCONJG(W(I,IPO))*AM
+ SUM2=SUM2+DCONJG(W(I,JPO))*AM
+ SUM3=SUM3+W(I,IPOPT)*AM
+ 12 SUM4=SUM4+W(I,JPOPT)*AM
+ SCSC(IPO)=CCCS*SUM
+ SCSCP(IPO)=CCCS*SUMP
+ ECSC(IPO)=-CCCS*DREAL(SUM1)
+ ECSCP(IPO)=-CCCS*SUM2
+ FSAC(IPO,IPO)=CSAM*SUM1
+ FSAC(JPO,IPO)=CSAM*SUM2
+ SAC(IPO,IPO)=CSAM*SUM3
+ 18 SAC(JPO,IPO)=CSAM*SUM4
+ I=0
+ DO 22 IPO1=1,2
+ DO 22 JPO1=1,2
+ CC=DCONJG(SAC(JPO1,IPO1))
+ DO 22 IPO2=1,2
+ DO 22 JPO2=1,2
+ I=I+1
+ 22 VINT(I)=SAC(JPO2,IPO2)*CC*CFSQ
+ RETURN
+ END
+ SUBROUTINE PCRSM0(VK,EXRI,INPOL)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+ COMPLEX*16 SUM1,SUM2,SUM3,SUM4,SUMPD,SUMS1,SUMS2,SUMS3,SUMS4,
+ 1CSAM,CC0,UIM
+ CC0=(0.0D0,0.0D0)
+ UIM=(0.0D0,1.0D0)
+ C0=0.0D0
+ EXDC=EXRI*EXRI
+ CCS=4.0D0*DACOS(C0)/(VK*VK)
+ CCCS=CCS/EXDC
+ CSAM=-(CCS/(EXRI*VK))*(0.0D0,0.5D0)
+ SUM2=CC0
+ SUM3=CC0
+ DO 14 I=1,NLEM
+ IE=I+NLEM
+ SUM2=SUM2+(AM0M(I,I)+AM0M(IE,IE))
+ 14 SUM3=SUM3+(AM0M(I,IE)+AM0M(IE,I))
+ SUMPI=C0
+ SUMPD=CC0
+ NLEMT=NLEM+NLEM
+ DO 16 I=1,NLEMT
+ DO 16 J=1,NLEM
+ JE=J+NLEM
+ SUMPI=SUMPI+DREAL
+ 1(DCONJG(AM0M(I,J))*AM0M(I,J)+DCONJG(AM0M(I,JE))*AM0M(I,JE))
+ 16 SUMPD=SUMPD+
+ 1(DCONJG(AM0M(I,J))*AM0M(I,JE)+DCONJG(AM0M(I,JE))*AM0M(I,J))
+ IF(INPOL.NE.0)GO TO 18
+ SUM1=SUM2
+ SUM4=SUM3*UIM
+ SUM3=-SUM4
+ SUMS1=SUMPI
+ SUMS2=SUMPI
+ SUMS3=SUMPD*UIM
+ SUMS4=-SUMS3
+ GO TO 20
+ 18 SUM1=SUM2+SUM3
+ SUM2=SUM2-SUM3
+ SUM3=CC0
+ SUM4=CC0
+ SUMS1=SUMPI-SUMPD
+ SUMS2=SUMPI+SUMPD
+ SUMS3=CC0
+ SUMS4=CC0
+ 20 ECSCM(1)=-CCCS*DREAL(SUM2)
+ ECSCM(2)=-CCCS*DREAL(SUM1)
+ ECSCPM(1)=-CCCS*SUM4
+ ECSCPM(2)=-CCCS*SUM3
+ FSACM(1,1)=CSAM*SUM2
+ FSACM(2,1)=CSAM*SUM4
+ FSACM(2,2)=CSAM*SUM1
+ FSACM(1,2)=CSAM*SUM3
+ SCSCM(1)=CCCS*DREAL(SUMS1)
+ SCSCM(2)=CCCS*DREAL(SUMS2)
+ SCSCPM(1)=CCCS*SUMS3
+ SCSCPM(2)=CCCS*SUMS4
+ RETURN
+ END
+ SUBROUTINE CRSM1(VK,EXRI)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+ COMMON/C6/RAC3J(17)
+CCC DIMENSION SVF(LE4PO,LE4PO,4),SVW(LE4PO,4,4),SVS(LE4PO,4)
+CCC LE4PO=LE*4+1
+ DIMENSION SVF(33,33,4),SVW(33,4,4),SVS(33,4)
+ COMPLEX*16 SVF,SVW,SVS
+ COMPLEX*16 CAM,CC0
+ CC0=(0.0D0,0.0D0)
+ C0=0.0D0
+ EXDC=EXRI*EXRI
+ CCS=1.0D0/(VK*VK)
+ PIGFSQ=6.4D+1*DACOS(C0)**2
+ CINT=CCS/(PIGFSQ*EXDC)
+ LETPO=LE+LE+1
+ DO 20 I=1,16
+ 20 VINTM(I)=CC0
+ DO 40 LPO=1,LETPO
+ L=LPO-1
+ LTPO=LPO+L
+ IMMN=LETPO-L
+ IMMX=LETPO+L
+ DO 22 IMF=IMMN,IMMX
+ DO 22 IMS=IMMN,IMMX
+ DO 22 IPO=1,4
+ 22 SVF(IMF,IMS,IPO)=CC0
+ DO 34 L1=1,LE
+ IL1=L1*(L1+1)
+ DO 34 L2=1,LE
+ IF((L.LT.IABS(L2-L1)).OR.(L.GT.(L2+L1)))GO TO 34
+ IL2=L2*(L2+1)
+ DO 24 IM=IMMN,IMMX
+ DO 24 IPA=1,4
+ SVS(IM,IPA)=CC0
+ DO 24 IPO=1,4
+ 24 SVW(IM,IPA,IPO)=CC0
+ DO 30 IM=IMMN,IMMX
+ M=IM-LETPO
+ CALL R3JMR(L,L1,L2,M)
+ M1MNMO=MAX0(-L1,-L2-M)-1
+ NM1=MIN0(L1,L2-M)-M1MNMO
+ DO 30 IM1=1,NM1
+ M1=-IM1-M1MNMO
+ ISN=1
+ IF(MOD(M1,2).NE.0)ISN=-1
+ CG3J=RAC3J(IM1)*ISN
+CCC CG=(LTPO**(1/2))*CG3J
+ ILM1=IL1+M1
+ ILM2=IL2+M1-M
+ IPA=0
+ DO 30 IPA1=1,2
+ I1=ILM1
+ IF(IPA1.EQ.2)I1=ILM1+NLEM
+ DO 30 IPA2=1,2
+ I2=ILM2
+ IF(IPA2.EQ.2)I2=ILM2+NLEM
+ IPA=IPA+1
+ SVS(IM,IPA)=SVS(IM,IPA)+AM0M(I1,I2)*CG3J
+ IPO=0
+ DO 30 IPO2=1,2
+ DO 30 IPO1=3,4
+ IPO=IPO+1
+ 30 SVW(IM,IPA,IPO)=SVW(IM,IPA,IPO)+W(I1,IPO1)*W(I2,IPO2)*CG3J
+ DO 32 IMF=IMMN,IMMX
+ DO 32 IMS=IMMN,IMMX
+ DO 32 IPO=1,4
+ DO 32 IPA=1,4
+ 32 SVF(IMF,IMS,IPO)=SVF(IMF,IMS,IPO)+SVW(IMF,IPA,IPO)*SVS(IMS,IPA)
+ 34 CONTINUE
+ DO 40 IMF=IMMN,IMMX
+ DO 40 IMS=IMMN,IMMX
+ I=0
+ DO 40 IPO1=1,4
+ CAM=DCONJG(SVF(IMF,IMS,IPO1))
+ DO 40 IPO2=1,4
+ I=I+1
+ 40 VINTM(I)=VINTM(I)+(SVF(IMF,IMS,IPO2)*CAM)*LTPO
+ DO 42 I=1,16
+ 42 VINTM(I)=VINTM(I)*CINT
+ RETURN
+ END
+ SUBROUTINE APS(ZPV,LI,NSPH,IOG,RMI,REI,SQK,GAPS)
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC DIMENSION ZPV(LM,3,2,2),IOG(NSPH),
+CCC 1RMI(LI,NSPH),REI(LI,NSPH),GAPS(NSPH)
+ DIMENSION ZPV(8,3,2,2),IOG(6),RMI(8,6),REI(8,6),GAPS(6)
+ COMPLEX*16 RMI,REI
+ COMPLEX*16 CC0,SUMM,SUME,SUEM,SUEE,SUM
+ CC0=(0.0D0,0.0D0)
+ C0=0.0D0
+ PIGH=DACOS(C0)
+ COFS=PIGH*2.0D0/SQK
+ DO 40 I=1,NSPH
+ IF(IOG(I).LT.I)GO TO 40
+ SUM=CC0
+ DO 30 L=1,LI
+ LTPO=L+L+1
+ DO 30 ILMP=1,3
+ IF(((L.EQ.1).AND.(ILMP.EQ.1)).OR.((L.EQ.LI).AND.(ILMP.EQ.3)))
+ 1GO TO 30
+ LMPML=ILMP-2
+ LMP=L+LMPML
+ COFL=LTPO*(LMP+LMP+1)
+ COFL=DSQRT(COFL)
+ SUMM=ZPV(L,ILMP,1,1)/(DCONJG(RMI(L,I))*RMI(LMP,I))
+ SUME=ZPV(L,ILMP,1,2)/(DCONJG(RMI(L,I))*REI(LMP,I))
+ SUEM=ZPV(L,ILMP,2,1)/(DCONJG(REI(L,I))*RMI(LMP,I))
+ SUEE=ZPV(L,ILMP,2,2)/(DCONJG(REI(L,I))*REI(LMP,I))
+ SUM=SUM+(CG1(LMPML,0,L,-1)*(SUMM-SUME-SUEM+SUEE)+
+ 1CG1(LMPML,0,L,1)*(SUMM+SUME+SUEM+SUEE))*COFL
+ 30 CONTINUE
+ GAPS(I)=DREAL(SUM)*COFS
+ 40 CONTINUE
+ RETURN
+ END
+ SUBROUTINE APC(ZPV,LE,AM0M,W,SQK,GAPR,GAPP)
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC DIMENSION ZPV(LM,3,2,2),AM0M(NLEMT,NLEMT),W(NLEMT,4)
+ DIMENSION ZPV(8,3,2,2),AM0M(160,160),W(160,4),
+ 1GAPR(3,2),GAPP(3,2)
+ COMPLEX*16 AM0M,W,GAPP
+CCC DIMENSION AC(NLEMT,2)
+ DIMENSION AC(160,2),GAP(3,2)
+ COMPLEX*16 AC,GAP
+ COMPLEX*16 CC0,UIM,UIMMP,
+ 1SUMM,SUME,SUEM,SUEE,SUMMP,SUMEP,SUEMP,SUEEP
+ CC0=(0.0D0,0.0D0)
+ UIM=(0.0D0,1.0D0)
+ COF=1.0D0/SQK
+ CIMU=COF/DSQRT(2.0D0)
+ NLEM=LE*(LE+2)
+ NLEMT=NLEM+NLEM
+ DO 45 J=1,NLEMT
+ AC(J,1)=CC0
+ AC(J,2)=CC0
+ DO 45 I=1,NLEMT
+ AC(J,1)=AC(J,1)+AM0M(J,I)*W(I,1)
+ AC(J,2)=AC(J,2)+AM0M(J,I)*W(I,2)
+ 45 CONTINUE
+ DO 90 IMU=1,3
+ MU=IMU-2
+ GAP(IMU,1)=CC0
+ GAP(IMU,2)=CC0
+ GAPP(IMU,1)=CC0
+ GAPP(IMU,2)=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)
+ JPO=2
+ DO 60 IPO=1,2
+ IF(IPO.EQ.2)JPO=1
+ SUMM=DCONJG(AC(I,IPO))*AC(IMP,IPO)
+ SUME=DCONJG(AC(I,IPO))*AC(IMPE,IPO)
+ SUEM=DCONJG(AC(IE,IPO))*AC(IMP,IPO)
+ SUEE=DCONJG(AC(IE,IPO))*AC(IMPE,IPO)
+ SUMMP=DCONJG(AC(I,JPO))*AC(IMP,IPO)
+ SUMEP=DCONJG(AC(I,JPO))*AC(IMPE,IPO)
+ SUEMP=DCONJG(AC(IE,JPO))*AC(IMP,IPO)
+ SUEEP=DCONJG(AC(IE,JPO))*AC(IMPE,IPO)
+ IF(LMPML.EQ.0)GO TO 55
+ SUMM=SUMM*UIMMP
+ SUME=SUME*UIMMP
+ SUEM=SUEM*UIMMP
+ SUEE=SUEE*UIMMP
+ SUMMP=SUMMP*UIMMP
+ SUMEP=SUMEP*UIMMP
+ SUEMP=SUEMP*UIMMP
+ SUEEP=SUEEP*UIMMP
+ 55 GAP(IMU,IPO)=GAP(IMU,IPO)+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))
+ GAPP(IMU,IPO)=GAPP(IMU,IPO)+CGC*
+ 1(SUMMP*ZPV(L,ILMP,1,1)+SUMEP*ZPV(L,ILMP,1,2)+
+ 2SUEMP*ZPV(L,ILMP,2,1)+SUEEP*ZPV(L,ILMP,2,2))
+ 60 CONTINUE
+ 70 CONTINUE
+ 80 CONTINUE
+ 90 CONTINUE
+ DO 95 IPO=1,2
+ SUME=GAP(1,IPO)*CIMU
+ SUEE=GAP(2,IPO)*COF
+ SUEM=GAP(3,IPO)*CIMU
+ GAPR(1,IPO)=DREAL(SUME-SUEM)
+ GAPR(2,IPO)=DREAL((SUME+SUEM)*UIM)
+ GAPR(3,IPO)=DREAL(SUEE)
+ SUMEP=GAPP(1,IPO)*CIMU
+ SUEEP=GAPP(2,IPO)*COF
+ SUEMP=GAPP(3,IPO)*CIMU
+ GAPP(1,IPO)=SUMEP-SUEMP
+ GAPP(2,IPO)=(SUMEP+SUEMP)*UIM
+ GAPP(3,IPO)=SUEEP
+ 95 CONTINUE
+ RETURN
+ END
+ SUBROUTINE APCRA(ZPV,LE,AM0M,INPOL,SQK,GAPRM,GAPPM)
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC DIMENSION ZPV(LM,3,2,2),AM0M(NLEMT,NLEMT)
+ DIMENSION ZPV(8,3,2,2),AM0M(160,160),GAPRM(3,2),GAPPM(3,2)
+ COMPLEX*16 AM0M,GAPPM
+CCC DIMENSION SVW(LE,3,2,2),SVS(LE,3,2,2)
+ DIMENSION SVW(8,3,2,2),SVS(8,3,2,2)
+ COMPLEX*16 SVS
+ COMPLEX*16 CC0,UIM,UIMTL,UIMTLS,FC,
+ 1CA11,CA12,CA21,CA22,A11,A12,A21,A22,SUM1,SUM2
+ CC0=(0.0D0,0.0D0)
+ UIM=(0.0D0,1.0D0)
+ NLEM=LE*(LE+2)
+ DO 28 L=1,LE
+ LPO=L+1
+ LTPO=LPO+L
+ FL=LTPO
+ FL=DSQRT(FL)
+ DO 26 ILMP=1,3
+ IF(((L.EQ.1).AND.(ILMP.EQ.1)).OR.((L.EQ.LE).AND.(ILMP.EQ.3)))
+ 1GO TO 26
+ LMPML=ILMP-2
+ LMP=L+LMPML
+ FLMP=LMP+LMP+1
+ FLMP=DSQRT(FLMP)
+ FLLMP=FLMP/FL
+ CGMMO=FLLMP*CG1(LMPML,0,L,1)
+ CGMPO=FLLMP*CG1(LMPML,0,L,-1)
+ IF(INPOL.NE.0)GO TO 22
+ CGS=CGMPO+CGMMO
+ CGD=CGMPO-CGMMO
+ SVW(L,ILMP,1,1)=CGS
+ SVW(L,ILMP,1,2)=CGD
+ SVW(L,ILMP,2,1)=CGD
+ SVW(L,ILMP,2,2)=CGS
+ GO TO 26
+ 22 SVW(L,ILMP,1,1)=CGMPO
+ SVW(L,ILMP,2,1)=CGMPO
+ SVW(L,ILMP,1,2)=-CGMMO
+ SVW(L,ILMP,2,2)=CGMMO
+ 26 CONTINUE
+ 28 CONTINUE
+ DO 30 L=1,LE
+ DO 30 ILMP=1,3
+ DO 30 IPA=1,2
+ DO 30 IPAMP=1,2
+ SVS(L,ILMP,IPA,IPAMP)=CC0
+ 30 CONTINUE
+ DO 58 L=1,LE
+ LPO=L+1
+ LTPO=L+LPO
+ IMM=L*LPO
+ DO 56 ILMP=1,3
+ IF(((L.EQ.1).AND.(ILMP.EQ.1)).OR.((L.EQ.LE).AND.(ILMP.EQ.3)))
+ 1GO TO 56
+ LMPML=ILMP-2
+ LMP=L+LMPML
+ IMPMM=LMP*(LMP+1)
+ UIMTL=LMPML*UIM
+ IF(LMPML.EQ.0)UIMTL=(1.0D0,0.0D0)
+ DO 54 IM=1,LTPO
+ M=IM-LPO
+ I=IMM+M
+ IE=I+NLEM
+ DO 52 IMU=1,3
+ MU=IMU-2
+ MMP=M-MU
+ IF(IABS(MMP).GT.LMP)GO TO 52
+ IMP=IMPMM+MMP
+ IMPE=IMP+NLEM
+ CGC=CG1(LMPML,-MU,L,-M)
+ DO 48 LS=1,LE
+ LSPO=LS+1
+ LSTPO=LS+LSPO
+ ISMM=LS*LSPO
+ DO 46 ILSMP=1,3
+ IF(((LS.EQ.1).AND.(ILSMP.EQ.1)).OR.((LS.EQ.LE).AND.(ILSMP.EQ.3)))
+ 1GO TO 46
+ LSMPML=ILSMP-2
+ LSMP=LS+LSMPML
+ ISMPMM=LSMP*(LSMP+1)
+ UIMTLS=-LSMPML*UIM
+ IF(LSMPML.EQ.0)UIMTLS=(1.0D0,0.0D0)
+ DO 44 IMS=1,LSTPO
+ MS=IMS-LSPO
+ MSMP=MS-MU
+ IF(IABS(MSMP).GT.LSMP)GO TO 44
+ IS=ISMM+MS
+ ISE=IS+NLEM
+ ISMP=ISMPMM+MSMP
+ ISMPE=ISMP+NLEM
+ CGCS=CG1(LSMPML,MU,LS,MS)
+ FC=(UIMTL*UIMTLS)*(CGC*CGCS)
+ CA11=DCONJG(AM0M(IS,I))
+ CA12=DCONJG(AM0M(IS,IE))
+ CA21=DCONJG(AM0M(ISE,I))
+ CA22=DCONJG(AM0M(ISE,IE))
+ A11=AM0M(ISMP,IMP)
+ A12=AM0M(ISMP,IMPE)
+ A21=AM0M(ISMPE,IMP)
+ A22=AM0M(ISMPE,IMPE)
+ Z11=ZPV(LS,ILSMP,1,1)
+ Z12=ZPV(LS,ILSMP,1,2)
+ Z21=ZPV(LS,ILSMP,2,1)
+ Z22=ZPV(LS,ILSMP,2,2)
+ SVS(L,ILMP,1,1)=SVS(L,ILMP,1,1)+
+ 1(CA11*A11*Z11+CA11*A21*Z12+CA21*A11*Z21+CA21*A21*Z22)*FC
+ SVS(L,ILMP,1,2)=SVS(L,ILMP,1,2)+
+ 1(CA11*A12*Z11+CA11*A22*Z12+CA21*A12*Z21+CA21*A22*Z22)*FC
+ SVS(L,ILMP,2,1)=SVS(L,ILMP,2,1)+
+ 1(CA12*A11*Z11+CA12*A21*Z12+CA22*A11*Z21+CA22*A21*Z22)*FC
+ SVS(L,ILMP,2,2)=SVS(L,ILMP,2,2)+
+ 1(CA12*A12*Z11+CA12*A22*Z12+CA22*A12*Z21+CA22*A22*Z22)*FC
+ 44 CONTINUE
+ 46 CONTINUE
+ 48 CONTINUE
+ 52 CONTINUE
+ 54 CONTINUE
+ 56 CONTINUE
+ 58 CONTINUE
+ SUM1=CC0
+ SUM2=CC0
+ DO 68 L=1,LE
+ LPO=L+1
+ LTPO=L+LPO
+ IMM=L*LPO
+ DO 66 ILMP=1,3
+ IF(((L.EQ.1).AND.(ILMP.EQ.1)).OR.((L.EQ.LE).AND.(ILMP.EQ.3)))
+ 1GO TO 66
+ IF(INPOL.NE.0)GO TO 62
+ SUM1=SUM1+
+ 1SVW(L,ILMP,1,1)*SVS(L,ILMP,1,1)+SVW(L,ILMP,2,1)*SVS(L,ILMP,1,2)+
+ 2SVW(L,ILMP,2,1)*SVS(L,ILMP,2,1)+SVW(L,ILMP,1,1)*SVS(L,ILMP,2,2)
+ SUM2=SUM2+
+ 1SVW(L,ILMP,1,2)*SVS(L,ILMP,1,1)+SVW(L,ILMP,2,2)*SVS(L,ILMP,1,2)+
+ 2SVW(L,ILMP,2,2)*SVS(L,ILMP,2,1)+SVW(L,ILMP,1,2)*SVS(L,ILMP,2,2)
+ GO TO 66
+ 62 SUM1=SUM1+
+ 1SVW(L,ILMP,2,1)*SVS(L,ILMP,1,1)+SVW(L,ILMP,1,1)*SVS(L,ILMP,1,2)+
+ 2SVW(L,ILMP,1,1)*SVS(L,ILMP,2,1)+SVW(L,ILMP,2,1)*SVS(L,ILMP,2,2)
+ SUM2=SUM2+
+ 1SVW(L,ILMP,2,2)*SVS(L,ILMP,1,1)+SVW(L,ILMP,1,2)*SVS(L,ILMP,1,2)+
+ 2SVW(L,ILMP,1,2)*SVS(L,ILMP,2,1)+SVW(L,ILMP,2,2)*SVS(L,ILMP,2,2)
+ 66 CONTINUE
+ 68 CONTINUE
+ C0=0.0D0
+ PIGH=DACOS(C0)
+ COFS=PIGH*2.0D0/SQK
+ GAPRM(1,1)=C0
+ GAPRM(1,2)=C0
+ GAPRM(2,1)=C0
+ GAPRM(2,2)=C0
+ GAPPM(1,1)=CC0
+ GAPPM(1,2)=CC0
+ GAPPM(2,1)=CC0
+ GAPPM(2,2)=CC0
+ IF(INPOL.NE.0)GO TO 72
+ SUM1=SUM1*COFS
+ SUM2=SUM2*COFS
+ GAPRM(3,1)=DREAL(SUM1)
+ GAPRM(3,2)=DREAL(SUM1)
+ GAPPM(3,1)=SUM2*UIM
+ GAPPM(3,2)=-GAPPM(3,1)
+ RETURN
+ 72 COFS=COFS*2.0D0
+ GAPRM(3,1)=DREAL(SUM1)*COFS
+ GAPRM(3,2)=DREAL(SUM2)*COFS
+ GAPPM(3,1)=CC0
+ GAPPM(3,2)=CC0
+ RETURN
+ END
+ SUBROUTINE THDPS(LM,ZPV)
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC DIMENSION ZPV(LM,3,2,2)
+ DIMENSION ZPV(8,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 UPVMP(THD,PHD,ICSPNV,COST,SINT,COSP,SINP,U,UP,UN)
+ IMPLICIT REAL*8(A-H,O-Z)
+ DIMENSION U(3),UP(3),UN(3)
+ C0=0.0D0
+ PIGH=DACOS(C0)
+ RDR=PIGH/9.0D+1
+ TH=THD*RDR
+ PH=PHD*RDR
+ COST=DCOS(TH)
+ SINT=DSIN(TH)
+ COSP=DCOS(PH)
+ SINP=DSIN(PH)
+ U(1)=COSP*SINT
+ U(2)=SINP*SINT
+ U(3)=COST
+ UP(1)=COSP*COST
+ UP(2)=SINP*COST
+ UP(3)=-SINT
+ UN(1)=-SINP
+ UN(2)=COSP
+ UN(3)=C0
+ IF(ICSPNV.EQ.0)RETURN
+ UP(1)=-UP(1)
+ UP(2)=-UP(2)
+ UP(3)=-UP(3)
+ UN(1)=-UN(1)
+ UN(2)=-UN(2)
+ RETURN
+ END
+ SUBROUTINE WMAMP
+ 1(IIS,COST,SINT,COSP,SINP,INPOL,LM,IDOT,NSPH,ARG,U,UP,UN)
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC DIMENSION YLM(NLMM+2),ARG(NSPEF)
+CCC NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE
+ DIMENSION YLM(82),U(3),UP(3),UN(3),ARG(1)
+ COMPLEX*16 YLM
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ NLMP=LM*(LM+2)+2
+ YLM(NLMP)=(0.0D0,0.0D0)
+ IF(IDOT.EQ.0)GO TO 65
+ IF(IDOT.EQ.1)GO TO 45
+ DO 40 N=1,NSPH
+ 40 ARG(N)=U(1)*RXX(N)+U(2)*RYY(N)+U(3)*RZZ(N)
+ GO TO 55
+ 45 DO 50 N=1,NSPH
+ 50 ARG(N)=COST*RZZ(N)
+ 55 IF(IIS.NE.2)GO TO 65
+ DO 60 N=1,NSPH
+ 60 ARG(N)=-ARG(N)
+ 65 CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM)
+ CALL PWMA(UP,UN,YLM,INPOL,LM,IIS)
+ RETURN
+ END
+ SUBROUTINE UPVSP(U,UPMP,UNMP,US,UPSMP,UNSMP,UP,UN,UPS,UNS,DUK,
+ 1ISQ,IBF,SCAND,CFMP,SFMP,CFSP,SFSP)
+ IMPLICIT REAL*8(A-H,O-Z)
+ DIMENSION U(3),UPMP(3),UNMP(3),US(3),UPSMP(3),UNSMP(3),
+ 1UP(3),UN(3),UPS(3),UNS(3),DUK(3)
+ C0=0.0D0
+ RDR=DACOS(C0)/9.0D+1
+ SML=1.0D-6
+ ISQ=0
+ SCAND=U(1)*US(1)+U(2)*US(2)+U(3)*US(3)
+ IF(DABS(SCAND-1.0D0).LT.SML)GO TO 10
+ IF(DABS(SCAND+1.0D0).LT.SML)GO TO 15
+ SCAND=DACOS(SCAND)/RDR
+ DUK(1)=U(1)-US(1)
+ DUK(2)=U(2)-US(2)
+ DUK(3)=U(3)-US(3)
+ IBF=0
+ GO TO 25
+ 10 SCAND=C0
+ DUK(1)=C0
+ DUK(2)=C0
+ DUK(3)=C0
+ IBF=-1
+ ISQ=-1
+ UPS(1)=UPSMP(1)
+ UPS(2)=UPSMP(2)
+ UPS(3)=UPSMP(3)
+ UNS(1)=UNSMP(1)
+ UNS(2)=UNSMP(2)
+ UNS(3)=UNSMP(3)
+ GO TO 20
+ 15 SCAND=180.0D0
+ DUK(1)=U(1)*2.0D0
+ DUK(2)=U(2)*2.0D0
+ DUK(3)=U(3)*2.0D0
+ IBF=1
+ UPS(1)=-UPSMP(1)
+ UPS(2)=-UPSMP(2)
+ UPS(3)=-UPSMP(3)
+ UNS(1)=-UNSMP(1)
+ UNS(2)=-UNSMP(2)
+ UNS(3)=-UNSMP(3)
+ 20 UP(1)=UPMP(1)
+ UP(2)=UPMP(2)
+ UP(3)=UPMP(3)
+ UN(1)=UNMP(1)
+ UN(2)=UNMP(2)
+ UN(3)=UNMP(3)
+ GO TO 85
+ 25 CALL ORUNVE(U,US,UN,-1,SML)
+ UNS(1)=UN(1)
+ UNS(2)=UN(2)
+ UNS(3)=UN(3)
+ CALL ORUNVE(UN,U,UP,1,SML)
+ CALL ORUNVE(UNS,US,UPS,1,SML)
+ 85 CFMP=UPMP(1)*UP(1)+UPMP(2)*UP(2)+UPMP(3)*UP(3)
+ SFMP=UNMP(1)*UP(1)+UNMP(2)*UP(2)+UNMP(3)*UP(3)
+ CFSP=UPS(1)*UPSMP(1)+UPS(2)*UPSMP(2)+UPS(3)*UPSMP(3)
+ SFSP=UNS(1)*UPSMP(1)+UNS(2)*UPSMP(2)+UNS(3)*UPSMP(3)
+ RETURN
+ END
+ SUBROUTINE WMASP(COST,SINT,COSP,SINP,COSTS,SINTS,COSPS,SINPS,
+ 1U,UP,UN,US,UPS,UNS,ISQ,IBF,INPOL,LM,IDOT,NSPH,ARGI,ARGS)
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC DIMENSION YLM(NLMM+2),ARGI(NSPEF),ARGS(NSPEF)
+CCC NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE
+ DIMENSION YLM(82),ARGI(1),ARGS(1),
+ 1U(3),UP(3),UN(3),US(3),UPS(3),UNS(3)
+ COMPLEX*16 YLM
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ NLMP=LM*(LM+2)+2
+ YLM(NLMP)=(0.0D0,0.0D0)
+ IF(IDOT.EQ.0)GO TO 75
+ IF(IDOT.EQ.1)GO TO 50
+ DO 40 N=1,NSPH
+ ARGI(N)=U(1)*RXX(N)+U(2)*RYY(N)+U(3)*RZZ(N)
+ IF(IBF.EQ.0)GO TO 35
+ ARGS(N)=ARGI(N)*IBF
+ GO TO 40
+ 35 ARGS(N)=-(US(1)*RXX(N)+US(2)*RYY(N)+US(3)*RZZ(N))
+ 40 CONTINUE
+ GO TO 75
+ 50 DO 60 N=1,NSPH
+ ARGI(N)=COST*RZZ(N)
+ IF(IBF.EQ.0)GO TO 55
+ ARGS(N)=ARGI(N)*IBF
+ GO TO 60
+ 55 ARGS(N)=-COSTS*RZZ(N)
+ 60 CONTINUE
+ 75 CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM)
+ CALL PWMA(UP,UN,YLM,INPOL,LM,ISQ)
+ IF(IBF.LT.0)RETURN
+ CALL SPHAR(COSTS,SINTS,COSPS,SINPS,LM,YLM)
+ CALL PWMA(UPS,UNS,YLM,INPOL,LM,2)
+ RETURN
+ END
+ SUBROUTINE STR(RCF)
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC DIMENSION RCF(NSPH,NSHL),YLM(MAX0(LITPOS,LMTPOS))
+ DIMENSION RCF(6,8),YLM(289)
+ COMPLEX*16 YLM
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ COMMON/C3/TFSAS,TSAS(2,2),GCS,SCS,ECS,ACS
+ COMPLEX*16 TFSAS,TSAS
+ COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+ COMMON/C6/RAC3J(17)
+ C0=0.0D0
+ PIGH=DACOS(C0)
+ PIG=PIGH*2.0D0
+ GCS=C0
+ DO 18 I=1,NSPH
+ IOGI=IOG(I)
+ IF(IOGI.LT.I)GO TO 18
+ GCSS=PIG*ROS(I)*ROS(I)
+ GCSV(I)=GCSS
+ NSH=NSHL(I)
+ DO 16 J=1,NSH
+ 16 RC(I,J)=RCF(I,J)*ROS(I)
+ 18 GCS=GCS+GCSV(IOGI)
+ LITPO=LI+LI+1
+ LITPOS=LITPO*LITPO
+ LMTPO=LI+LE+1
+ LMTPOS=LMTPO*LMTPO
+ NLIM=LI*(LI+2)
+ NLEM=LE*(LE+2)
+ LM=MAX0(LI,LE)
+ LMPO=LM+1
+ I=0
+ DO 28 L1PO=1,LMPO
+ L1=L1PO-1
+ DO 28 L2=1,LM
+ CALL R3J000(L1,L2)
+ IND3J(L1PO,L2)=I
+ LMNPO=IABS(L2-L1)+1
+ LMXPO=L2+L1+1
+ IL=0
+ DO 28 LPO=LMNPO,LMXPO,2
+ I=I+1
+ IL=IL+1
+ 28 V3J0(I)=RAC3J(IL)
+ NSPHMO=NSPH-1
+ LIT=LI+LI
+ IVY=0
+ DO 40 NF=1,NSPHMO
+ NFPO=NF+1
+ DO 40 NS=NFPO,NSPH
+ RX=RXX(NF)-RXX(NS)
+ RY=RYY(NF)-RYY(NS)
+ RZ=RZZ(NF)-RZZ(NS)
+ CALL POLAR(RX,RY,RZ,RR,CRTH,SRTH,CRPH,SRPH)
+ CALL SPHAR(CRTH,SRTH,CRPH,SRPH,LIT,YLM)
+ DO 38 IV=1,LITPOS
+ 38 VYHJ(IV+IVY)=DCONJG(YLM(IV))
+ 40 IVY=IVY+LITPOS
+ LMT=LI+LE
+ IVY=0
+ DO 50 NF=1,NSPH
+ RX=RXX(NF)
+ RY=RYY(NF)
+ RZ=RZZ(NF)
+ IF(RX.EQ.C0.AND.RY.EQ.C0.AND.RZ.EQ.C0)GO TO 50
+ CALL POLAR(RX,RY,RZ,RR,CRTH,SRTH,CRPH,SRPH)
+ CALL SPHAR(CRTH,SRTH,CRPH,SRPH,LMT,YLM)
+ DO 48 IV=1,LMTPOS
+ 48 VYJ0(IV+IVY)=DCONJG(YLM(IV))
+ 50 IVY=IVY+LMTPOS
+ RETURN
+ END
+ SUBROUTINE HJV(EXRI,VK,JER,LCALC,ARG)
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC DIMENSION RFJ(MAX0(LIT,LMT)),RFN(LITPO)
+ DIMENSION RFJ(17),RFN(17)
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+ COMPLEX*16 ARG,UIM
+ JER=0
+ C0=0.0D0
+ UIM=(0.0D0,1.0D0)
+ NSPHMO=NSPH-1
+ LIT=LI+LI
+ IVHB=0
+ DO 40 NF=1,NSPHMO
+ NFPO=NF+1
+ DO 40 NS=NFPO,NSPH
+ RX=RXX(NF)-RXX(NS)
+ RY=RYY(NF)-RYY(NS)
+ RZ=RZZ(NF)-RZZ(NS)
+ RR=DSQRT(RX*RX+RY*RY+RZ*RZ)
+ RARG=RR*VK*EXRI
+ ARG=RARG
+ CALL RBF(LIT,RARG,LCALC,RFJ)
+ IF(LCALC.GE.LIT)GO TO 35
+ JER=1
+ RETURN
+ 35 CALL RNF(LIT,RARG,LCALC,RFN)
+ IF(LCALC.GE.LIT)GO TO 37
+ JER=2
+ RETURN
+ 37 DO 38 LPO=1,LITPO
+ 38 VH(LPO+IVHB)=RFJ(LPO)+UIM*RFN(LPO)
+ 40 IVHB=IVHB+LITPO
+ LMT=LI+LE
+ IVHB=0
+ DO 50 NF=1,NSPH
+ RX=RXX(NF)
+ RY=RYY(NF)
+ RZ=RZZ(NF)
+ IF(RX.EQ.C0.AND.RY.EQ.C0.AND.RZ.EQ.C0)GO TO 50
+ RR=DSQRT(RX*RX+RY*RY+RZ*RZ)
+ RARG=RR*VK*EXRI
+ CALL RBF(LMT,RARG,LCALC,RFJ)
+ IF(LCALC.GE.LMT)GO TO 45
+ JER=3
+ ARG=RARG
+ RETURN
+ 45 DO 47 LPO=1,LMTPO
+ 47 VJ0(LPO+IVHB)=RFJ(LPO)
+ 50 IVHB=IVHB+LMTPO
+ RETURN
+ END
+ SUBROUTINE PWMA(UP,UN,YLM,INPOL,LW,ISQ)
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC DIMENSION YLM(NLWM+2)
+ DIMENSION YLM(82),UP(3),UN(3)
+ COMPLEX*16 YLM
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ COMPLEX*16 CP1,CM1,C1,CP2,CM2,C2,UIM,CL
+ PI4=DACOS(0.0D0)*8.0D0
+ IS=ISQ
+ IF(ISQ.EQ.-1)IS=0
+ ISPO=IS+1
+ ISPT=IS+2
+ NLWM=LW*(LW+2)
+ NLWMT=NLWM+NLWM
+ SQRTWI=1.0D0/DSQRT(2.0D0)
+ 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,ISPO)=DCONJG(CP1*CP*YLM(K+2)+CM1*CM*YLM(K)+CZ1*CZ*YLM(K+1))*CL
+ 20 W(K,ISPT)=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,ISPO)=UIM*W(K,ISPT)
+ 30 W(I,ISPT)=-UIM*W(K,ISPO)
+ 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
+ W(I,ISPO)=-C2
+ W(K,ISPT)=C1
+ 40 W(I,ISPT)=C1
+ 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))
+ RETURN
+ END
+ SUBROUTINE ORUNVE(U1,U2,U3,IORTH,TORTH)
+ IMPLICIT REAL*8(A-H,O-Z)
+ DIMENSION U1(3),U2(3),U3(3)
+ IF(IORTH.GT.0)GO TO 10
+ CP=U1(1)*U2(1)+U1(2)*U2(2)+U1(3)*U2(3)
+ IF((IORTH.LT.0).AND.(DABS(CP).LT.TORTH))GO TO 10
+ FN=1.0D0/DSQRT(1.0D0-CP*CP)
+ U3(1)=(U1(2)*U2(3)-U1(3)*U2(2))*FN
+ U3(2)=(U1(3)*U2(1)-U1(1)*U2(3))*FN
+ U3(3)=(U1(1)*U2(2)-U1(2)*U2(1))*FN
+ RETURN
+ 10 U3(1)=U1(2)*U2(3)-U1(3)*U2(2)
+ U3(2)=U1(3)*U2(1)-U1(1)*U2(3)
+ U3(3)=U1(1)*U2(2)-U1(2)*U2(1)
+ RETURN
+ END
+ SUBROUTINE CMS(AM)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+CCC DIMENSION AM(MXNDM,MXNDM)
+CCC NDIT.LE.MXNDM
+ DIMENSION AM(960,960)
+ COMPLEX*16 AM
+ COMPLEX*16 DM,DE,CGH,CGK,GHIT,CC0
+ CC0=(0.0D0,0.0D0)
+ NDI=NSPH*NLIM
+ NBL=0
+ NSPHMO=NSPH-1
+ DO 26 N1=1,NSPHMO
+ IN1=(N1-1)*NLIM
+ N1PO=N1+1
+ DO 26 N2=N1PO,NSPH
+ IN2=(N2-1)*NLIM
+ NBL=NBL+1
+ DO 24 L1=1,LI
+ L1PO=L1+1
+ IL1=L1PO*L1
+ L1TPO=L1PO+L1
+ DO 24 IM1=1,L1TPO
+ M1=IM1-L1PO
+ ILM1=IL1+M1
+ ILM1E=ILM1+NDI
+ I1=IN1+ILM1
+ I1E=IN1+ILM1E
+ J1=IN2+ILM1
+ J1E=IN2+ILM1E
+ DO 24 L2=1,LI
+ L2PO=L2+1
+ IL2=L2PO*L2
+ L2TPO=L2PO+L2
+ ISH=(-1)**(L2+L1)
+ ISK=-ISH
+ DO 24 IM2=1,L2TPO
+ M2=IM2-L2PO
+ ILM2=IL2+M2
+ ILM2E=ILM2+NDI
+ I2=IN2+ILM2
+ I2E=IN2+ILM2E
+ J2=IN1+ILM2
+ J2E=IN1+ILM2E
+ CGH=GHIT(0,0,NBL,L1,M1,L2,M2)
+ CGK=GHIT(0,1,NBL,L1,M1,L2,M2)
+ AM(I1,I2)=CGH
+ AM(I1,I2E)=CGK
+ AM(I1E,I2)=CGK
+ AM(I1E,I2E)=CGH
+ AM(J1,J2)=CGH*ISH
+ AM(J1,J2E)=CGK*ISK
+ AM(J1E,J2)=CGK*ISK
+ AM(J1E,J2E)=CGH*ISH
+ 24 CONTINUE
+ 26 CONTINUE
+ DO 30 N1=1,NSPH
+ IN1=(N1-1)*NLIM
+ DO 30 L1=1,LI
+ DM=RMI(L1,N1)
+ DE=REI(L1,N1)
+ L1PO=L1+1
+ IL1=L1PO*L1
+ L1TPO=L1PO+L1
+ DO 30 IM1=1,L1TPO
+ M1=IM1-L1PO
+ ILM1=IL1+M1
+ I1=IN1+ILM1
+ I1E=I1+NDI
+ DO 28 ILM2=1,NLIM
+ I2=IN1+ILM2
+ I2E=I2+NDI
+ AM(I1,I2)=CC0
+ AM(I1,I2E)=CC0
+ AM(I1E,I2)=CC0
+ 28 AM(I1E,I2E)=CC0
+ AM(I1,I1)=DM
+ 30 AM(I1E,I1E)=DE
+ RETURN
+ END
+ SUBROUTINE ZTM(AM)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+ COMMON/C9/GIS(480,80),GLS(480,80),SAM(960,160)
+ COMPLEX*16 GIS,GLS,SAM
+CCC DIMENSION AM(MXNDM,MXNDM)
+CCC NDIT.LE.MXNDM
+ DIMENSION AM(960,960)
+ COMPLEX*16 AM
+ COMPLEX*16 GIE,GLE,GHIT,A1,A2,A3,A4,
+ 1SUM1,SUM2,SUM3,SUM4,CC0
+ CC0=(0.0D0,0.0D0)
+ NDI=NSPH*NLIM
+ I2=0
+ DO 15 N2=1,NSPH
+ DO 15 L2=1,LI
+ L2TPO=L2+L2+1
+ M2=-L2-1
+ DO 15 IM2=1,L2TPO
+ M2=M2+1
+ I2=I2+1
+ I3=0
+ DO 15 L3=1,LE
+ L3TPO=L3+L3+1
+ M3=-L3-1
+ DO 15 IM3=1,L3TPO
+ M3=M3+1
+ I3=I3+1
+ GIS(I2,I3)=GHIT(2,0,N2,L2,M2,L3,M3)
+ 15 GLS(I2,I3)=GHIT(2,1,N2,L2,M2,L3,M3)
+ DO 25 I1=1,NDI
+ I1E=I1+NDI
+ DO 25 I3=1,NLEM
+ I3E=I3+NLEM
+ SUM1=CC0
+ SUM2=CC0
+ SUM3=CC0
+ SUM4=CC0
+ DO 20 I2=1,NDI
+ I2E=I2+NDI
+ GIE=GIS(I2,I3)
+ GLE=GLS(I2,I3)
+ A1=AM(I1,I2)
+ A2=AM(I1,I2E)
+ A3=AM(I1E,I2)
+ A4=AM(I1E,I2E)
+ SUM1=SUM1+A1*GIE+A2*GLE
+ SUM2=SUM2+A1*GLE+A2*GIE
+ SUM3=SUM3+A3*GIE+A4*GLE
+ 20 SUM4=SUM4+A3*GLE+A4*GIE
+ SAM(I1,I3)=SUM1
+ SAM(I1,I3E)=SUM2
+ SAM(I1E,I3)=SUM3
+ 25 SAM(I1E,I3E)=SUM4
+ DO 35 I1=1,NDI
+ DO 35 I0=1,NLEM
+ GIS(I1,I0)=DCONJG(GIS(I1,I0))
+ 35 GLS(I1,I0)=DCONJG(GLS(I1,I0))
+ NLEMT=NLEM+NLEM
+ DO 45 I0=1,NLEM
+ I0E=I0+NLEM
+ DO 45 I3=1,NLEMT
+ SUM1=CC0
+ SUM2=CC0
+ DO 40 I1=1,NDI
+ I1E=I1+NDI
+ A1=SAM(I1,I3)
+ A2=SAM(I1E,I3)
+ GIE=GIS(I1,I0)
+ GLE=GLS(I1,I0)
+ SUM1=SUM1+A1*GIE+A2*GLE
+ 40 SUM2=SUM2+A1*GLE+A2*GIE
+ AM0M(I0,I3)=-SUM1
+ 45 AM0M(I0E,I3)=-SUM2
+ RETURN
+ END
+ COMPLEX*16 FUNCTION GHIT(IHI,IPAMO,NBL,L1,M1,L2,M2)
+CCC NBL INDIVIDUATES TRANSFER VECTOR GOING FROM N2 TO N1;
+CCC IHI=0 FOR HANKEL, IHI=1 FOR BESSEL, IHI=2 FOR BESSEL FROM ORIGIN;
+CCC DEPENDING ON IHI, IPAMO=0 GIVES H OR I, IPAMO=1 GIVES K OR L
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+ COMMON/C6/RAC3J(17)
+ COMPLEX*16 CSUM,CC0,UIM,CFUN
+ CC0=(0.0D0,0.0D0)
+ GHIT=CC0
+ IF(IHI.NE.2)GO TO 10
+ C0=0.0D0
+ IF(RXX(NBL).NE.C0.OR.RYY(NBL).NE.C0.OR.RZZ(NBL).NE.C0)GO TO 10
+ IF(IPAMO.NE.0)RETURN
+ IF(L1.EQ.L2.AND.M1.EQ.M2)GHIT=(1.0D0,0.0D0)
+ RETURN
+ 10 L1MP=L1-IPAMO
+ L1PO=L1+1
+ M1MM2=M1-M2
+ M1MM2M=IABS(M1MM2)+1
+ LMINPO=IABS(L2-L1MP)+1
+ LMAXPO=L2+L1MP+1
+ I3J0IN=IND3J(L1MP+1,L2)
+ ILIN=-1
+ IF((M1MM2M.GT.LMINPO).AND.(MOD(M1MM2M-LMINPO,2).NE.0))ILIN=0
+ ISN=1
+ IF(MOD(M1,2).NE.0)ISN=-ISN
+ IF(MOD(LMINPO,2).NE.0)GO TO 12
+ ISN=-ISN
+ IF(L2.GT.L1MP)ISN=-ISN
+ 12 NBLMO=NBL-1
+ IF(IHI.EQ.2)GO TO 50
+ NBHJ=NBLMO*LITPO
+ NBY=NBLMO*LITPOS
+ IF(IHI.EQ.1)GO TO 30
+ DO 24 JM=1,3
+ CSUM=CC0
+ MU=JM-2
+ MUPM1=MU+M1
+ MUPM2=MU+M2
+ IF(MUPM1.LT.-L1MP.OR.MUPM1.GT.L1MP.OR.MUPM2.LT.-L2.OR.MUPM2.GT.L2)
+ 1GO TO 24
+ JSN=-ISN
+ IF(MU.EQ.0)JSN=ISN
+ CR=CGEV(IPAMO,MU,L1,M1)*CGEV(0,MU,L2,M2)
+ I3J0=I3J0IN
+ IF((MUPM1.NE.0).OR.(MUPM2.NE.0))GO TO 16
+ DO 14 LT=LMINPO,LMAXPO,2
+ I3J0=I3J0+1
+ L3=LT-1
+ NY=L3*L3+LT
+ AORS=L3+LT
+ F3J=(V3J0(I3J0)*V3J0(I3J0)*DSQRT(AORS))*JSN
+ CFUN=(VH(NBHJ+LT)*VYHJ(NBY+NY))*F3J
+ CSUM=CSUM+CFUN
+ 14 JSN=-JSN
+ GO TO 22
+ 16 CALL R3JJR(L1MP,L2,-MUPM1,MUPM2)
+ IL=ILIN
+ DO 20 LT=LMINPO,LMAXPO,2
+ I3J0=I3J0+1
+ IF(M1MM2M.GT.LT)GO TO 20
+ IL=IL+2
+ L3=LT-1
+ NY=L3*L3+LT+M1MM2
+ AORS=L3+LT
+ F3J=(RAC3J(IL)*V3J0(I3J0)*DSQRT(AORS))*JSN
+ CFUN=(VH(NBHJ+LT)*VYHJ(NBY+NY))*F3J
+ CSUM=CSUM+CFUN
+ 20 JSN=-JSN
+ 22 CSUM=CSUM*CR
+ GHIT=GHIT+CSUM
+ 24 CONTINUE
+ GO TO 70
+ 30 DO 44 JM=1,3
+ CSUM=CC0
+ MU=JM-2
+ MUPM1=MU+M1
+ MUPM2=MU+M2
+ IF(MUPM1.LT.-L1MP.OR.MUPM1.GT.L1MP.OR.MUPM2.LT.-L2.OR.MUPM2.GT.L2)
+ 1GO TO 44
+ JSN=-ISN
+ IF(MU.EQ.0)JSN=ISN
+ CR=CGEV(IPAMO,MU,L1,M1)*CGEV(0,MU,L2,M2)
+ I3J0=I3J0IN
+ IF((MUPM1.NE.0).OR.(MUPM2.NE.0))GO TO 36
+ DO 34 LT=LMINPO,LMAXPO,2
+ I3J0=I3J0+1
+ L3=LT-1
+ NY=L3*L3+LT
+ AORS=L3+LT
+ F3J=(V3J0(I3J0)*V3J0(I3J0)*DSQRT(AORS))*JSN
+ CFUN=(VJ(NBHJ+LT)*VYHJ(NBY+NY))*F3J
+ CSUM=CSUM+CFUN
+ 34 JSN=-JSN
+ GO TO 42
+ 36 CALL R3JJR(L1MP,L2,-MUPM1,MUPM2)
+ IL=ILIN
+ DO 40 LT=LMINPO,LMAXPO,2
+ I3J0=I3J0+1
+ IF(M1MM2M.GT.LT)GO TO 40
+ IL=IL+2
+ L3=LT-1
+ NY=L3*L3+LT+M1MM2
+ AORS=L3+LT
+ F3J=(RAC3J(IL)*V3J0(I3J0)*DSQRT(AORS))*JSN
+ CFUN=(VJ(NBHJ+LT)*VYHJ(NBY+NY))*F3J
+ CSUM=CSUM+CFUN
+ 40 JSN=-JSN
+ 42 CSUM=CSUM*CR
+ GHIT=GHIT+CSUM
+ 44 CONTINUE
+ GO TO 70
+ 50 NBHJ=NBLMO*LMTPO
+ NBY=NBLMO*LMTPOS
+ DO 64 JM=1,3
+ CSUM=CC0
+ MU=JM-2
+ MUPM1=MU+M1
+ MUPM2=MU+M2
+ IF(MUPM1.LT.-L1MP.OR.MUPM1.GT.L1MP.OR.MUPM2.LT.-L2.OR.MUPM2.GT.L2)
+ 1GO TO 64
+ JSN=-ISN
+ IF(MU.EQ.0)JSN=ISN
+ CR=CGEV(IPAMO,MU,L1,M1)*CGEV(0,MU,L2,M2)
+ I3J0=I3J0IN
+ IF((MUPM1.NE.0).OR.(MUPM2.NE.0))GO TO 56
+ DO 54 LT=LMINPO,LMAXPO,2
+ I3J0=I3J0+1
+ L3=LT-1
+ NY=L3*L3+LT
+ AORS=L3+LT
+ F3J=(V3J0(I3J0)*V3J0(I3J0)*DSQRT(AORS))*JSN
+ CFUN=(VJ0(NBHJ+LT)*VYJ0(NBY+NY))*F3J
+ CSUM=CSUM+CFUN
+ 54 JSN=-JSN
+ GO TO 62
+ 56 CALL R3JJR(L1MP,L2,-MUPM1,MUPM2)
+ IL=ILIN
+ DO 60 LT=LMINPO,LMAXPO,2
+ I3J0=I3J0+1
+ IF(M1MM2M.GT.LT)GO TO 60
+ IL=IL+2
+ L3=LT-1
+ NY=L3*L3+LT+M1MM2
+ AORS=L3+LT
+ F3J=(RAC3J(IL)*V3J0(I3J0)*DSQRT(AORS))*JSN
+ CFUN=(VJ0(NBHJ+LT)*VYJ0(NBY+NY))*F3J
+ CSUM=CSUM+CFUN
+ 60 JSN=-JSN
+ 62 CSUM=CSUM*CR
+ GHIT=GHIT+CSUM
+ 64 CONTINUE
+ 70 PI4=DACOS(0.0D0)*8.0D0
+ IF(IPAMO.EQ.1)GO TO 90
+ CR=PI4*(L1+L1PO)*(L2+L2+1)
+ CR=DSQRT(CR)
+ GHIT=GHIT*CR
+ RETURN
+ 90 UIM=(0.0D0,1.0D0)
+ CR=L1PO
+ CR=(PI4*(L1+L1MP)*(L1+L1PO)*(L2+L2+1))/CR
+ CR=DSQRT(CR)
+ GHIT=GHIT*(CR*UIM)
+ RETURN
+ END
+ REAL*8 FUNCTION CGEV(IPAMO,MU,L,M)
+CCC CGEV(IPAMO,MU,L,M)=CLGO(1,L-IPAMO,L;-MU,MU+M,M)
+ REAL*8 XD,XN
+ IF(IPAMO.NE.0)GO TO 30
+ IF((M.NE.0).OR.(MU.NE.0))GO TO 10
+ CGEV=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)
+ CGEV=DSQRT(XN/XD)
+ RETURN
+ 15 XN=(L-M)*(L+M+1)
+ CGEV=-DSQRT(XN/XD)
+ RETURN
+ 20 XD=(L+1)*L
+ XN=-M
+ CGEV=XN/DSQRT(XD)
+ RETURN
+ 30 XD=(L*2-1)*L*2
+ IF(MU)35,40,45
+ 35 XN=(L-1+M)*(L+M)
+ GO TO 50
+ 40 XN=(L-M)*(L+M)*2
+ GO TO 50
+ 45 XN=(L-1-M)*(L-M)
+ 50 CGEV=DSQRT(XN/XD)
+ RETURN
+ END
+ SUBROUTINE R3J000(J2,J3)
+CCC 3j(J,J2,J3;0,0,0)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C6/RAC3J(17)
+ ONE=1.0D0
+ JMX=J3+J2
+ IF(JMX.GT.0)GO TO 10
+ RAC3J(1)=ONE
+ RETURN
+ 10 JMN=IABS(J3-J2)
+ NJMO=(JMX-JMN)/2
+ JF=JMX+JMX+1
+ ISN=1
+ IF(MOD(JMN,2).NE.0)ISN=-1
+ IF(NJMO.GT.0)GO TO 15
+ SJ=JF
+ CNR=(ONE/DSQRT(SJ))*ISN
+ RAC3J(1)=CNR
+ RETURN
+ 15 SJR=JF
+ JMXPOS=(JMX+1)*(JMX+1)
+ JMNS=JMN*JMN
+ J1MO=JMX-1
+ J1S=(J1MO+1)*(J1MO+1)
+ CJ=(JMXPOS-J1S)*(J1S-JMNS)
+ CJ=DSQRT(CJ)
+ J1MOS=J1MO*J1MO
+ CJMO=(JMXPOS-J1MOS)*(J1MOS-JMNS)
+ CJMO=DSQRT(CJMO)
+ IF(NJMO.GT.1)GO TO 20
+ RAC3J(1)=-CJ/CJMO
+ SJ=SJR+(RAC3J(1)*RAC3J(1))*(JF-4)
+ CNR=(ONE/DSQRT(SJ))*ISN
+ RAC3J(2)=CNR
+ RAC3J(1)=RAC3J(1)*CNR
+ RETURN
+ 20 NJ=NJMO+1
+ NMAT=(NJ+1)/2
+ RAC3J(NJ)=ONE
+ RAC3J(NJMO)=-CJ/CJMO
+ IF(NMAT.EQ.NJMO)GO TO 50
+ NBR=NJMO-NMAT
+ DO 45 IBR=1,NBR
+ IRR=NJ-IBR
+ JF=JF-4
+ J1MO=J1MO-2
+ J1S=(J1MO+1)*(J1MO+1)
+ CJ=(JMXPOS-J1S)*(J1S-JMNS)
+ CJ=DSQRT(CJ)
+ J1MOS=J1MO*J1MO
+ CJMO=(JMXPOS-J1MOS)*(J1MOS-JMNS)
+ CJMO=DSQRT(CJMO)
+ RAC3J(IRR-1)=RAC3J(IRR)*(-CJ/CJMO)
+ 45 SJR=SJR+(RAC3J(IRR)*RAC3J(IRR))*JF
+ 50 RACMAT=RAC3J(NMAT)
+ SJR=SJR+(RACMAT*RACMAT)*(JF-4)
+ RAC3J(1)=ONE
+ JF=JMN+JMN+1
+ SJL=JF
+ J1PT=JMN+2
+ J1POS=(J1PT-1)*(J1PT-1)
+ CJPO=(JMXPOS-J1POS)*(J1POS-JMNS)
+ CJPO=DSQRT(CJPO)
+ J1PTS=J1PT*J1PT
+ CJPT=(JMXPOS-J1PTS)*(J1PTS-JMNS)
+ CJPT=DSQRT(CJPT)
+ RAC3J(2)=-CJPO/CJPT
+ NMATMO=NMAT-1
+ IF(NMATMO.LT.2)GO TO 75
+ DO 70 IRL=2,NMATMO
+ JF=JF+4
+ J1PT=J1PT+2
+ J1POS=(J1PT-1)*(J1PT-1)
+ CJPO=(JMXPOS-J1POS)*(J1POS-JMNS)
+ CJPO=DSQRT(CJPO)
+ J1PTS=J1PT*J1PT
+ CJPT=(JMXPOS-J1PTS)*(J1PTS-JMNS)
+ CJPT=DSQRT(CJPT)
+ RAC3J(IRL+1)=RAC3J(IRL)*(-CJPO/CJPT)
+ 70 SJL=SJL+(RAC3J(IRL)*RAC3J(IRL))*JF
+ 75 RATRAC=RACMAT/RAC3J(NMAT)
+ RATS=RATRAC*RATRAC
+ SJ=SJR+SJL*RATS
+ RAC3J(NMAT)=RACMAT
+ CNR=(ONE/DSQRT(SJ))*ISN
+ DO 80 IRR=NMAT,NJ
+ 80 RAC3J(IRR)=RAC3J(IRR)*CNR
+ CNL=CNR*RATRAC
+ DO 85 IRL=1,NMATMO
+ 85 RAC3J(IRL)=RAC3J(IRL)*CNL
+ RETURN
+ END
+ SUBROUTINE R3JJR(J2,J3,M2,M3)
+CCC 3j(J,J2,J3;-M2-M3,M2,M3)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C6/RAC3J(17)
+ ONE=1.0D0
+ JMX=J3+J2
+ JDF=J3-J2
+ M1=-M2-M3
+ JMN=MAX0(IABS(JDF),IABS(M1))
+ NJMO=JMX-JMN
+ JF=JMX+JMX+1
+ ISN=1
+ IF(MOD(JDF+M1,2).NE.0)ISN=-1
+ IF(NJMO.GT.0)GO TO 15
+ SJ=JF
+ CNR=(ONE/DSQRT(SJ))*ISN
+ RAC3J(1)=CNR
+ RETURN
+ 15 SJT=ONE
+ SJR=JF
+ JSMPOS=(JMX+1)*(JMX+1)
+ JDFS=JDF*JDF
+ M1S=M1*M1
+ MDF=M3-M2
+ IDJC=M1*(J3*(J3+1)-J2*(J2+1))
+ J1=JMX
+ J1S=J1*J1
+ J1PO=J1+1
+ CCJ=(J1S-JDFS)*(J1S-M1S)
+ CJ=CCJ*(JSMPOS-J1S)
+ CJ=DSQRT(CJ)
+ DJ=JF*(J1*J1PO*MDF+IDJC)
+ IF(NJMO.GT.1)GO TO 20
+ RAC3J(1)=-DJ/(CJ*J1PO)
+ SJ=SJR+(RAC3J(1)*RAC3J(1))*(JF-2)
+ CNR=(ONE/DSQRT(SJ))*ISN
+ RAC3J(2)=CNR
+ RAC3J(1)=RAC3J(1)*CNR
+ RETURN
+ 20 NJ=NJMO+1
+ NMAT=(NJ+1)/2
+ RAC3J(NJ)=ONE
+ RAC3J(NJMO)=-DJ/(CJ*J1PO)
+ IF(NMAT.EQ.NJMO)GO TO 50
+ NBR=NJMO-NMAT
+ DO 45 IBR=1,NBR
+ IRR=NJ-IBR
+ JF=JF-2
+ J1=J1-1
+ J1S=J1*J1
+ J1PO=J1+1
+ CJP=CJ
+ CCJ=(J1S-JDFS)*(J1S-M1S)
+ CJ=CCJ*(JSMPOS-J1S)
+ CJ=DSQRT(CJ)
+ SJT=RAC3J(IRR)*RAC3J(IRR)
+ DJ=JF*(J1*J1PO*MDF+IDJC)
+ RAC3J(IRR-1)=-(RAC3J(IRR)*DJ+RAC3J(IRR+1)*(CJP*J1))/(CJ*J1PO)
+ 45 SJR=SJR+SJT*JF
+ 50 OSJT=SJT
+ SJT=RAC3J(NMAT)*RAC3J(NMAT)
+ IF(SJT.LT.OSJT)GO TO 55
+ SJR=SJR+SJT*(JF-2)
+ GO TO 60
+ 55 NMAT=NMAT+1
+ 60 RACMAT=RAC3J(NMAT)
+ RAC3J(1)=ONE
+ JF=JMN+JMN+1
+ SJL=JF
+ J1=JMN
+ IF(J1.EQ.0)GO TO 62
+ J1PO=J1+1
+ J1POS=J1PO*J1PO
+ CCJP=(J1POS-JDFS)*(J1POS-M1S)
+ CJP=CCJP*(JSMPOS-J1POS)
+ CJP=DSQRT(CJP)
+ DJ=JF*(J1*J1PO*MDF+IDJC)
+ RAC3J(2)=-DJ/(CJP*J1)
+ GO TO 63
+ 62 CJP=JSMPOS-1
+ CJP=DSQRT(CJP)
+ DJ=MDF
+ RAC3J(2)=-DJ/CJP
+ 63 NMATMO=NMAT-1
+ IF(NMATMO.LT.2)GO TO 75
+ DO 70 IRL=2,NMATMO
+ JF=JF+2
+ J1=J1+1
+ J1PO=J1+1
+ J1POS=J1PO*J1PO
+ CJ=CJP
+ CCJP=(J1POS-JDFS)*(J1POS-M1S)
+ CJP=CCJP*(JSMPOS-J1POS)
+ CJP=DSQRT(CJP)
+ SJT=RAC3J(IRL)*RAC3J(IRL)
+ DJ=JF*(J1*J1PO*MDF+IDJC)
+ RAC3J(IRL+1)=-(RAC3J(IRL)*DJ+RAC3J(IRL-1)*(CJ*J1PO))/(CJP*J1)
+ 70 SJL=SJL+SJT*JF
+ 75 RATRAC=RACMAT/RAC3J(NMAT)
+ RATS=RATRAC*RATRAC
+ SJ=SJR+SJL*RATS
+ RAC3J(NMAT)=RACMAT
+ CNR=(ONE/DSQRT(SJ))*ISN
+ DO 80 IRR=NMAT,NJ
+ 80 RAC3J(IRR)=RAC3J(IRR)*CNR
+ CNL=CNR*RATRAC
+ DO 85 IRL=1,NMATMO
+ 85 RAC3J(IRL)=RAC3J(IRL)*CNL
+ RETURN
+ END
+ SUBROUTINE R3JMR(J1,J2,J3,M1)
+CCC 3j(J1,J2,J3;M1,M,-M1-M)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C6/RAC3J(17)
+ ONE=1.0D0
+ MMX=MIN0(J2,J3-M1)
+ MMN=MAX0(-J2,-J3-M1)
+ NMMO=MMX-MMN
+ J1PO=J1+1
+ J1TPO=J1PO+J1
+ ISN=1
+ IF(MOD(J2-J3-M1,2).NE.0)ISN=-1
+ IF(NMMO.GT.0)GO TO 15
+ SJ=J1TPO
+ CNR=(ONE/DSQRT(SJ))*ISN
+ RAC3J(1)=CNR
+ RETURN
+ 15 J1S=J1*J1PO
+ J2PO=J2+1
+ J2S=J2*J2PO
+ J3PO=J3+1
+ J3S=J3*J3PO
+ ID=J1S-J2S-J3S
+ M2=MMX
+ M3=M1+M2
+ CM=(J2PO-M2)*(J2+M2)*(J3PO-M3)*(J3+M3)
+ CM=DSQRT(CM)
+ DM=ID+M2*M3*2
+ IF(NMMO.GT.1)GO TO 20
+ RAC3J(1)=DM/CM
+ SJ=(ONE+RAC3J(1)*RAC3J(1))*J1TPO
+ CNR=(ONE/DSQRT(SJ))*ISN
+ RAC3J(2)=CNR
+ RAC3J(1)=RAC3J(1)*CNR
+ RETURN
+ 20 NM=NMMO+1
+ NMAT=(NM+1)/2
+ RAC3J(NM)=ONE
+ RAC3J(NMMO)=DM/CM
+ SJT=ONE
+ SJR=ONE
+ IF(NMAT.EQ.NMMO)GO TO 50
+ NBR=NMMO-NMAT
+ DO 45 IBR=1,NBR
+ IRR=NM-IBR
+ M2=M2-1
+ M3=M1+M2
+ CMP=CM
+ CM=(J2PO-M2)*(J2+M2)*(J3PO-M3)*(J3+M3)
+ CM=DSQRT(CM)
+ SJT=RAC3J(IRR)*RAC3J(IRR)
+ DM=ID+M2*M3*2
+ RAC3J(IRR-1)=(RAC3J(IRR)*DM-RAC3J(IRR+1)*CMP)/CM
+ 45 SJR=SJR+SJT
+ 50 OSJT=SJT
+ SJT=RAC3J(NMAT)*RAC3J(NMAT)
+ IF(SJT.LT.OSJT)GO TO 55
+ SJR=SJR+SJT
+ GO TO 60
+ 55 NMAT=NMAT+1
+ 60 RACMAT=RAC3J(NMAT)
+ RAC3J(1)=ONE
+ M2=MMN
+ M3=M1+M2
+ CMP=(J2-M2)*(J2PO+M2)*(J3-M3)*(J3PO+M3)
+ CMP=DSQRT(CMP)
+ DM=ID+M2*M3*2
+ RAC3J(2)=DM/CMP
+ SJL=ONE
+ NMATMO=NMAT-1
+ IF(NMATMO.LT.2)GO TO 75
+ DO 70 IRL=2,NMATMO
+ M2=M2+1
+ M3=M1+M2
+ CM=CMP
+ CMP=(J2-M2)*(J2PO+M2)*(J3-M3)*(J3PO+M3)
+ CMP=DSQRT(CMP)
+ SJT=RAC3J(IRL)*RAC3J(IRL)
+ DM=ID+M2*M3*2
+ RAC3J(IRL+1)=(RAC3J(IRL)*DM-RAC3J(IRL-1)*CM)/CMP
+ 70 SJL=SJL+SJT
+ 75 RATRAC=RACMAT/RAC3J(NMAT)
+ RATS=RATRAC*RATRAC
+ SJ=(SJR+SJL*RATS)*J1TPO
+ RAC3J(NMAT)=RACMAT
+ CNR=(ONE/DSQRT(SJ))*ISN
+ DO 80 IRR=NMAT,NM
+ 80 RAC3J(IRR)=RAC3J(IRR)*CNR
+ CNL=CNR*RATRAC
+ DO 85 IRL=1,NMATMO
+ 85 RAC3J(IRL)=RAC3J(IRL)*CNL
+ RETURN
+ END
+ SUBROUTINE DME(LI,I,NPNT,NPNTTS,VK,EXDC,EXRI,JER,LCALC,ARG)
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC DIMENSION CFJ(LIPT),RFJ(LIPT),RFN(LIPT),
+CCC 1FBI(LIPT),FB(LIPT),FN(LIPT),
+CCC 2RMF(LI),DRMF(LI),REF(LI),DREF(LI)
+ DIMENSION CFJ(10),RFJ(10),RFN(10),
+ 1FBI(10),FB(10),FN(10),RMF(8),DRMF(8),REF(8),DREF(8)
+ COMPLEX*16 CFJ,FBI,FB,FN,RMF,DRMF,REF,DREF
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ COMMON/C2/RIS(999),DLRI(999),DC0(5),VKT(6),VSZ(6)
+ COMPLEX*16 RIS,DLRI,DC0,VKT
+ COMPLEX*16 DFBI,DFB,DFN,CCNA,CCNB,CCNC,CCND,
+ 1Y1,DY1,Y2,DY2,ARIN,ARG,CRI,UIM
+ JER=0
+ UIM=(0.0D0,1.0D0)
+ NSTP=NPNT-1
+ NSTPTS=NPNTTS-1
+ LIPO=LI+1
+ LIPT=LI+2
+ SZ=VK*ROS(I)
+ VSZ(I)=SZ
+ VKR1=VK*RC(I,1)
+ NSH=NSHL(I)
+ VKT(I)=CDSQRT(DC0(1))
+ ARG=VKR1*VKT(I)
+ ARIN=ARG
+ IF(DIMAG(ARG).EQ.0.0D0)GO TO 26
+ CALL CBF(LIPO,ARG,LCALC,CFJ)
+ IF(LCALC.GE.LIPO)GO TO 22
+ JER=5
+ RETURN
+ 22 DO 24 J=1,LIPT
+ 24 FBI(J)=CFJ(J)
+ GO TO 32
+ 26 RARG=DREAL(ARG)
+ CALL RBF(LIPO,RARG,LCALC,RFJ)
+ IF(LCALC.GE.LIPO)GO TO 28
+ JER=5
+ RETURN
+ 28 DO 30 J=1,LIPT
+ 30 FBI(J)=RFJ(J)
+ 32 AREX=SZ*EXRI
+ ARG=AREX
+ RARG=AREX
+ CALL RBF(LIPO,RARG,LCALC,RFJ)
+ IF(LCALC.GE.LIPO)GO TO 40
+ JER=7
+ RETURN
+ 40 CALL RNF(LIPO,RARG,LCALC,RFN)
+ IF(LCALC.GE.LIPO)GO TO 42
+ JER=8
+ RETURN
+ 42 DO 43 J=1,LIPT
+ FB(J)=RFJ(J)
+ 43 FN(J)=RFN(J)
+ IF(NSH.GT.1)GO TO 65
+ CRI=DC0(1)/EXDC
+ DO 60 L=1,LI
+ LPO=L+1
+ LTPO=LPO+L
+ LPT=LPO+1
+ DFBI=(L*FBI(L)-LPO*FBI(LPT))*ARIN+FBI(LPO)*LTPO
+ DFB=(L*FB(L)-LPO*FB(LPT))*AREX+FB(LPO)*LTPO
+ DFN=(L*FN(L)-LPO*FN(LPT))*AREX+FN(LPO)*LTPO
+ CCNA=FBI(LPO)*DFN
+ CCNB=FN(LPO)*DFBI
+ CCNC=FBI(LPO)*DFB
+ CCND=FB(LPO)*DFBI
+ RMI(L,I)=1.0D0+UIM*(CCNA-CCNB)/(CCNC-CCND)
+ 60 REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND)
+ RETURN
+ 65 DO 80 L=1,LI
+ LPO=L+1
+ LTPO=LPO+L
+ LPT=LPO+1
+ DLTPO=LTPO
+cccccccccccccccccccccc
+ Y1=FBI(LPO)
+ DY1=(L*FBI(L)-LPO*FBI(LPT))*VKT(I)/DLTPO
+cccccccccccccccccccccc
+ Y2=Y1
+ DY2=DY1
+ IC=1
+ DO 76 NS=2,NSH
+ NSMO=NS-1
+ VKR=VK*RC(I,NSMO)
+ IF(MOD(NS,2).EQ.0)GO TO 70
+ IC=IC+1
+ STEP=NSTP
+ STEP=VK*(RC(I,NS)-RC(I,NSMO))/STEP
+ ARG=DC0(IC)
+ CALL RKC(NSTP,STEP,ARG,VKR,LPO,Y1,Y2,DY1,DY2)
+ GO TO 76
+ 70 CALL DIEL(NSTPTS,NSMO,I,IC,VK)
+ STEPTS=NSTPTS
+ STEPTS=VK*(RC(I,NS)-RC(I,NSMO))/STEPTS
+ CALL RKT(NSTPTS,STEPTS,VKR,LPO,Y1,Y2,DY1,DY2)
+ 76 CONTINUE
+ RMF(L)=Y1*SZ
+ DRMF(L)=DY1*SZ+Y1
+ REF(L)=Y2*SZ
+ 80 DREF(L)=DY2*SZ+Y2
+ CRI=(1.0D0,0.0D0)
+ IF(MOD(NSH,2).NE.0)CRI=DC0(IC)/EXDC
+ DO 90 L=1,LI
+ LPO=L+1
+ LTPO=LPO+L
+ LPT=LPO+1
+ DFB=(L*FB(L)-LPO*FB(LPT))*AREX+FB(LPO)*LTPO
+ DFN=(L*FN(L)-LPO*FN(LPT))*AREX+FN(LPO)*LTPO
+ CCNA=RMF(L)*DFN
+ CCNB=DRMF(L)*FN(LPO)*SZ*LTPO
+ CCNC=RMF(L)*DFB
+ CCND=DRMF(L)*FB(LPO)*SZ*LTPO
+ RMI(L,I)=1.0D0+UIM*(CCNA-CCNB)/(CCNC-CCND)
+ 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)
+ RETURN
+ END
+ SUBROUTINE RKT(NPNTMO,STEP,X,LPO,Y1,Y2,DY1,DY2)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C2/RIS(999),DLRI(999),DC0(5),VKT(6),VSZ(6)
+ COMPLEX*16 RIS,DLRI,DC0,VKT
+ COMPLEX*16 Y1,DY1,Y2,DY2,CY1,CDY1,C11,
+ 1CY23,CDY23,YC2,C12,C13,CY4,CDY4,YY,C14,C21,C22,C23,C24
+ HSTEP=.5D0*STEP
+ CL=LPO*(LPO-1)
+ DO 60 IPNT=1,NPNTMO
+ JPNT=IPNT+IPNT-1
+ CY1=CL/(X*X)-RIS(JPNT)
+ CDY1=-2.0D0/X
+ C11=(CY1*Y1+CDY1*DY1)*STEP
+ XH=X+HSTEP
+ JPNTPO=JPNT+1
+ CY23=CL/(XH*XH)-RIS(JPNTPO)
+ CDY23=-2.0D0/XH
+ YC2=Y1+DY1*HSTEP
+ C12=(CY23*YC2+CDY23*(DY1+.5D0*C11))*STEP
+ C13=(CY23*(YC2+.25D0*C11*STEP)+CDY23*(DY1+.5D0*C12))*STEP
+ XN=X+STEP
+ JPNTPT=JPNT+2
+ CY4=CL/(XN*XN)-RIS(JPNTPT)
+ CDY4=-2.0D0/XN
+ YY=Y1+DY1*STEP
+ C14=(CY4*(YY+.5D0*C12*STEP)+CDY4*(DY1+C13))*STEP
+ Y1=YY+(C11+C12+C13)*STEP/6.0D0
+ DY1=DY1+(.5D0*C11+C12+C13+.5D0*C14)/3.0D0
+ CY1=CY1-CDY1*DLRI(JPNT)
+ CDY1=CDY1+2.0D0*DLRI(JPNT)
+ C21=(CY1*Y2+CDY1*DY2)*STEP
+ CY23=CY23-CDY23*DLRI(JPNTPO)
+ CDY23=CDY23+2.0D0*DLRI(JPNTPO)
+ YC2=Y2+DY2*HSTEP
+ C22=(CY23*YC2+CDY23*(DY2+.5D0*C21))*STEP
+ C23=(CY23*(YC2+.25D0*C21*STEP)+CDY23*(DY2+.5D0*C22))*STEP
+ CY4=CY4-CDY4*DLRI(JPNTPT)
+ CDY4=CDY4+2.0D0*DLRI(JPNTPT)
+ YY=Y2+DY2*STEP
+ C24=(CY4*(YC2+.5D0*C22*STEP)+CDY4*(DY2+C23))*STEP
+ Y2=YY+(C21+C22+C23)*STEP/6.0D0
+ DY2=DY2+(.5D0*C21+C22+C23+.5D0*C24)/3.0D0
+ 60 X=XN
+ RETURN
+ END
+ SUBROUTINE RKC(NPNTMO,STEP,DCC,X,LPO,Y1,Y2,DY1,DY2)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMPLEX*16 Y1,DY1,Y2,DY2,CY1,CDY1,C11,
+ 1CY23,YC2,C12,C13,CY4,YY,C14,C21,C22,C23,C24,DCC
+ HSTEP=.5D0*STEP
+ CL=LPO*(LPO-1)
+ DO 60 IPNT=1,NPNTMO
+ CY1=CL/(X*X)-DCC
+ CDY1=-2.0D0/X
+ C11=(CY1*Y1+CDY1*DY1)*STEP
+ XH=X+HSTEP
+ CY23=CL/(XH*XH)-DCC
+ CDY23=-2.0D0/XH
+ YC2=Y1+DY1*HSTEP
+ C12=(CY23*YC2+CDY23*(DY1+.5D0*C11))*STEP
+ C13=(CY23*(YC2+.25D0*C11*STEP)+CDY23*(DY1+.5D0*C12))*STEP
+ XN=X+STEP
+ CY4=CL/(XN*XN)-DCC
+ CDY4=-2.0D0/XN
+ YY=Y1+DY1*STEP
+ C14=(CY4*(YY+.5D0*C12*STEP)+CDY4*(DY1+C13))*STEP
+ Y1=YY+(C11+C12+C13)*STEP/6.0D0
+ DY1=DY1+(.5D0*C11+C12+C13+.5D0*C14)/3.0D0
+ C21=(CY1*Y2+CDY1*DY2)*STEP
+ YC2=Y2+DY2*HSTEP
+ C22=(CY23*YC2+CDY23*(DY2+.5D0*C21))*STEP
+ C23=(CY23*(YC2+.25D0*C21*STEP)+CDY23*(DY2+.5D0*C22))*STEP
+ YY=Y2+DY2*STEP
+ C24=(CY4*(YC2+.5D0*C22*STEP)+CDY4*(DY2+C23))*STEP
+ Y2=YY+(C21+C22+C23)*STEP/6.0D0
+ DY2=DY2+(.5D0*C21+C22+C23+.5D0*C24)/3.0D0
+ 60 X=XN
+ RETURN
+ END
+ SUBROUTINE DIEL(NPNTMO,NS,I,IC,VK)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+ 1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+ 2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+ 3FSAC(2,2),SAC(2,2),FSACM(2,2),
+ 4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+ 5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+ 6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+ 7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+ 8IND3J(9,8),IOG(6),NSHL(6)
+ COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+ 1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+ COMMON/C2/RIS(999),DLRI(999),DC0(5),VKT(6),VSZ(6)
+ COMPLEX*16 RIS,DLRI,DC0,VKT
+ COMPLEX*16 DEL
+ DIF=RC(I,NS+1)-RC(I,NS)
+ HSTEP=NPNTMO
+ HSTEP=.5D0*DIF/HSTEP
+ RR=RC(I,NS)
+ DEL=DC0(IC+1)-DC0(IC)
+ KPNT=NPNTMO+NPNTMO
+ RIS(KPNT+1)=DC0(IC+1)
+ DLRI(KPNT+1)=(0.0D0,0.0D0)
+ DO 90 NP=1,KPNT
+ FF=(RR-RC(I,NS))/DIF
+ RIS(NP)=DEL*FF*FF*(-2.0D0*FF+3.0D0)+DC0(IC)
+ DLRI(NP)=3.0D0*DEL*FF*(1.0D0-FF)/(DIF*VK*RIS(NP))
+ 90 RR=RR+HSTEP
+ RETURN
+ END
+ SUBROUTINE RBF(N,X,NM,SJ)
+C
+C FROM SPHJ OF LIBRARY specfun
+C
+C =======================================================
+C Purpose: Compute spherical Bessel functions j
+C Input : X --- Argument of j
+C N --- Order of j ( N = 0,1,2,... )
+C Output: SJ(N+1) --- j
+C NM --- Highest order computed
+C Routines called:
+C MSTA1 and MSTA2 for computing the starting
+C point for backward recurrence
+C =======================================================
+C
+ IMPLICIT REAL*8 (A-H,O-Z)
+ DIMENSION SJ(N+1)
+ A0=DABS(X)
+ NM=N
+ IF(A0.LT.1.0D-60)THEN
+ DO K=2,N+1
+ SJ(K)=0.0D0
+ ENDDO
+ SJ(1)=1.0D0
+ RETURN
+ ENDIF
+ SJ(1)=DSIN(X)/X
+ IF(N.EQ.0)RETURN
+ SJ(2)=(SJ(1)-DCOS(X))/X
+ IF(N.EQ.1)RETURN
+ SA=SJ(1)
+ SB=SJ(2)
+ M=MSTA1(A0,200)
+ IF(M.LT.N)THEN
+ NM=M
+ ELSE
+ M=MSTA2(A0,N,15)
+ ENDIF
+ F0=0.0D0
+ F1=1.0D0-100
+ DO K=M,0,-1
+ F=(2.0D0*K+3.0D0)*F1/X-F0
+ IF(K.LE.NM)SJ(K+1)=F
+ F0=F1
+ F1=F
+ ENDDO
+ IF(DABS(SA).GT.DABS(SB))CS=SA/F
+ IF(DABS(SA).LE.DABS(SB))CS=SB/F0
+ DO K=1,NM+1
+ SJ(K)=CS*SJ(K)
+ ENDDO
+ END
+ SUBROUTINE CBF(N,Z,NM,CSJ)
+C
+C FROM CSPHJY OF LIBRARY specfun
+C
+C ==========================================================
+C Purpose: Compute spherical Bessel functions j
+C Input : Z --- Complex argument of j
+C N --- Order of j ( N = 0,1,2,... )
+C Output: CSJ(N+1) --- j
+C NM --- Highest order computed
+C Routines called:
+C MSTA1 and MSTA2 for computing the starting
+C point for backward recurrence
+C ==========================================================
+C
+ IMPLICIT COMPLEX*16 (C,Z)
+ REAL*8 A0
+ DIMENSION CSJ(N+1)
+ A0=CDABS(Z)
+ NM=N
+ IF(A0.LT.1.0D-60)THEN
+ DO K=2,N+1
+ CSJ(K)=0.0D0
+ ENDDO
+ CSJ(1)=(1.0D0,0.0D0)
+ RETURN
+ ENDIF
+ CSJ(1)=CDSIN(Z)/Z
+ IF(N.EQ.0)RETURN
+ CSJ(2)=(CSJ(1)-CDCOS(Z))/Z
+ IF(N.EQ.1)RETURN
+ CSA=CSJ(1)
+ CSB=CSJ(2)
+ M=MSTA1(A0,200)
+ IF(M.LT.N)THEN
+ NM=M
+ ELSE
+ M=MSTA2(A0,N,15)
+ ENDIF
+ CF0=0.0D0
+ CF1=1.0D0-100
+ DO K=M,0,-1
+ CF=(2.0D0*K+3.0D0)*CF1/Z-CF0
+ IF(K.LE.NM)CSJ(K+1)=CF
+ CF0=CF1
+ CF1=CF
+ ENDDO
+ IF(CDABS(CSA).GT.CDABS(CSB))CS=CSA/CF
+ IF(CDABS(CSA).LE.CDABS(CSB))CS=CSB/CF0
+ DO K=1,NM+1
+ CSJ(K)=CS*CSJ(K)
+ ENDDO
+ END
+ INTEGER FUNCTION MSTA1(X,MP)
+C
+C ===================================================
+C Purpose: Determine the starting point for backward
+C recurrence such that the magnitude of all
+C J or j at that point is about 10**(-MP)
+C Input : X --- Abs of argument of J or j
+C MP --- Value of magnitude
+C Output: MSTA1 --- Starting point
+C ===================================================
+C
+ IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+ A0=DABS(X)
+ N0=INT(1.1*A0)+1
+ F0=ENVJ(N0,A0)-MP
+ N1=N0+5
+ F1=ENVJ(N1,A0)-MP
+ DO 10 IT=1,20
+ NN=N1-IDINT((N1-N0)/(1.0D0-F0/F1))
+ F=ENVJ(NN,A0)-MP
+ IF(ABS(NN-N1).LT.1)GO TO 20
+ N0=N1
+ F0=F1
+ N1=NN
+ 10 F1=F
+ 20 MSTA1=NN
+ RETURN
+ END
+ INTEGER FUNCTION MSTA2(X,N,MP)
+C
+C ===================================================
+C Purpose: Determine the starting point for backward
+C recurrence such that all J or j
+C have MP significant digits
+C Input : X --- Abs of argument of J or j
+C N --- Order of J or j
+C MP --- Significant digit
+C Output: MSTA2 --- Starting point
+C ===================================================
+C
+ IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+ A0=DABS(X)
+ HMP=0.5D0*MP
+ EJN=ENVJ(N,A0)
+ IF(EJN.LE.HMP)THEN
+ OBJ=MP
+C
+ N0=INT(1.1*A0)+1
+C
+ ELSE
+ OBJ=HMP+EJN
+ N0=N
+ ENDIF
+ F0=ENVJ(N0,A0)-OBJ
+ N1=N0+5
+ F1=ENVJ(N1,A0)-OBJ
+ DO 10 IT=1,20
+ NN=N1-IDINT((N1-N0)/(1.0D0-F0/F1))
+ F=ENVJ(NN,A0)-OBJ
+ IF(ABS(NN-N1).LT.1)GO TO 20
+ N0=N1
+ F0=F1
+ N1=NN
+ 10 F1=F
+ 20 MSTA2=NN+10
+ RETURN
+ END
+ REAL*8 FUNCTION ENVJ(N,X)
+ DOUBLE PRECISION X,XN
+C
+ IF(N.EQ.0)THEN
+ XN=1.0D-100
+ ENVJ=0.5D0*DLOG10(6.28D0*XN)-XN*DLOG10(1.36D0*X/XN)
+ ELSE
+C
+ ENVJ=0.5D0*DLOG10(6.28D0*N)-N*DLOG10(1.36D0*X/N)
+C
+ ENDIF
+C
+ RETURN
+ END
+ SUBROUTINE RNF(N,X,NM,SY)
+C
+C FROM SPHJY OF LIBRARY specfun
+C
+C
+C ======================================================
+C Purpose: Compute spherical Bessel functions y
+C Input : X --- Argument of y ( X > 0 )
+C N --- Order of y ( N = 0,1,2,... )
+C Output: SY(N+1) --- y
+C NM --- Highest order computed
+C ======================================================
+C
+ IMPLICIT REAL*8 (A-H,O-Z)
+ DIMENSION SY(N+1)
+ IF(X.LT.1.0D-60)THEN
+ DO K=1,N+1
+ SY(K)=-1.0D+300
+ ENDDO
+ RETURN
+ ENDIF
+ SY(1)=-DCOS(X)/X
+ IF(N.EQ.0)RETURN
+ SY(2)=(SY(1)-DSIN(X))/X
+ IF(N.EQ.1)RETURN
+ F0=SY(1)
+ F1=SY(2)
+ DO K=2,N
+ F=(2.0D0*K-1.0D0)*F1/X-F0
+ SY(K+1)=F
+ IF(DABS(F).GE.1.0D+300)GO TO 20
+ F0=F1
+ F1=F
+ ENDDO
+ RETURN
+ 20 NM=K
+ RETURN
+ END
+ SUBROUTINE POLAR(X,Y,Z,R,CTH,STH,CPH,SPH)
+ IMPLICIT REAL*8(A-H,O-Z)
+ LOGICAL ONX,ONY,ONZ
+ 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 15
+ IF(ONY)GO TO 20
+ RHOS=X*X+Y*Y
+ RHO=DSQRT(RHOS)
+ CPH=X/RHO
+ SPH=Y/RHO
+ GO TO 25
+ 10 CPH=ONE
+ SPH=C0
+ GO TO 25
+ 15 RHOS=X*X
+ RHO=DABS(X)
+ CPH=DSIGN(ONE,X)
+ SPH=C0
+ GO TO 25
+ 20 RHOS=Y*Y
+ RHO=DABS(Y)
+ CPH=C0
+ SPH=DSIGN(ONE,Y)
+ 25 IF(Z.NE.C0)GO TO 35
+ IF(ONZ)GO TO 30
+ R=RHO
+ CTH=C0
+ STH=ONE
+ RETURN
+ 30 R=C0
+ CTH=ONE
+ STH=C0
+ RETURN
+ 35 IF(ONZ)GO TO 40
+ R=DSQRT(RHOS+Z*Z)
+ CTH=Z/R
+ STH=RHO/R
+ RETURN
+ 40 R=DABS(Z)
+ CTH=DSIGN(ONE,Z)
+ STH=C0
+ 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)
+CCC LL=MAX0(LIT,LMT)
+ DIMENSION SINRMP(16),COSRMP(16),PLEGN(153),YLM(289)
+ 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 LUCIN(AM,NDDMST,N,IER)
+C
+C NDDMST FIRST DIMENSION OF AM AS DECLARED IN DIMENSION
+C STATEMENT.
+C N NUMBER OF ROWS IN AM.
+C IER IS REPLACED BY 1 FOR SINGULARITY.
+C
+C
+C CALLS COMPLEX INNER PRODUCT FUNCTION CDTP.
+C
+ IMPLICIT REAL*8 (A-H,O-Z)
+CCC DIMENSION AM(NDDMST,NDDMST),V(NDDMST)
+ DIMENSION AM(960,960),V(960)
+ COMPLEX*16 AM
+ COMPLEX*16 CTEMP,CDTP
+ IER=0
+ NMINUS=N-1
+C
+C 1.2 DEAL WITH THE SPECIAL CASE N=1.
+C
+C1200 IF(N.NE.1)GO TO 1300
+C CTEMP=AM(1,1)
+C IF(CTEMP.EQ.(0.0D0,0.0D0))GO TO 5200
+C AM(1,1)=1.0D0/CTEMP
+C RETURN
+C
+C 1.3 SET V(I)=1/R(I)**2, WHERE R(I) IS THE LENGTH OF ROW I.
+C
+C1300 DO 1309 I=1,N
+ DO 1309 I=1,N
+ SUM=0.0D0
+ DO 1319 J=1,N
+ SUM=SUM+DREAL(AM(I,J))**2+DIMAG(AM(I,J))**2
+ 1319 CONTINUE
+ V(I)=1.0D0/SUM
+ 1309 CONTINUE
+C
+C 2. REPLACE AM BY TRIANGULAR MATRICES (L,U) WHERE AM=L*U.
+C REPLACE L(I,I) BY 1/L(I,I), READY FOR SECTION 4.
+C (ROW INTERCHANGES TAKE PLACE, AND THE INDICES OF THE PIVOTAL ROWS
+C ARE PLACED IN V.)
+C
+ DO 2019 K=1,N
+ KPLUS=K+1
+ KMINUS=K-1
+ L=K
+ PSQMAX=0.0D0
+ DO 2029 I=K,N
+ CTEMP=-CDTP(-AM(I,K),AM,I,1,K,KMINUS,NDDMST)
+ AM(I,K)=CTEMP
+ PSQ=V(I)*(DREAL(CTEMP)**2+DIMAG(CTEMP)**2)
+ IF(PSQ.LE.PSQMAX)GO TO 2029
+ PSQMAX=PSQ
+ L=I
+ 2029 CONTINUE
+ IF(L.EQ.K)GO TO 2011
+ DO 2049 J=1,N
+ CTEMP=AM(K,J)
+ AM(K,J)=AM(L,J)
+ AM(L,J)=CTEMP
+ 2049 CONTINUE
+ V(L)=V(K)
+ 2011 V(K)=L
+ IF(PSQMAX.EQ.0.0D0)GO TO 5200
+ CTEMP=1.0D0/AM(K,K)
+ AM(K,K)=CTEMP
+ IF(KPLUS.GT.N)GO TO 2019
+ DO 2059 J=KPLUS,N
+ AM(K,J)=-CTEMP*CDTP(-AM(K,J),AM,K,1,J,KMINUS,NDDMST)
+ 2059 CONTINUE
+ 2019 CONTINUE
+C
+C 4. REPLACE AM BY ITS INVERSE AMINV.
+C
+C
+C 4.1 REPLACE L AND U BY THEIR INVERSES LINV AND UINV.
+C
+ DO 4109 K=1,NMINUS
+ KPLUS=K+1
+ DO 4119 I=KPLUS,N
+ AM(I,K)=-AM(I,I)*CDTP((0.0D0,0.0D0),AM,I,K,K,I-K,
+ 1NDDMST)
+ AM(K,I)=-CDTP(AM(K,I),AM,K,KPLUS,I,I-K-1,NDDMST)
+ 4119 CONTINUE
+ 4109 CONTINUE
+C
+C 4.2 FORM AMINV=UINV*LINV.
+C
+ DO 4209 K=1,N
+ DO 4219 I=1,N
+ IF(I.GE.K)GO TO 4212
+ AM(I,K)=CDTP((0.0D0,0.0D0),AM,I,K,K,N-K+1,NDDMST)
+ GO TO 4219
+ 4212 AM(I,K)=CDTP(AM(I,K),AM,I,I+1,K,N-I,NDDMST)
+ 4219 CONTINUE
+ 4209 CONTINUE
+C
+C 4.3 INTERCHANGE COLUMNS OF AMINV AS SPECIFIED BY V, BUT IN REVERSE
+C ORDER.
+C
+ DO 4309 L=1,N
+ K=N-L+1
+ KCOL=IDINT(V(K))
+ IF(KCOL.EQ.K)GO TO 4309
+ DO 4319 I=1,N
+ CTEMP=AM(I,K)
+ AM(I,K)=AM(I,KCOL)
+ AM(I,KCOL)=CTEMP
+ 4319 CONTINUE
+ 4309 CONTINUE
+ RETURN
+ 5200 IER=1
+ RETURN
+ END
+ COMPLEX*16 FUNCTION CDTP(Z,AV,I,JF,K,NJ,ISTEP)
+CCC DIMENSION AV(NDDMST*NDDMST)
+ DIMENSION AV(921600)
+ COMPLEX*16 AV,Z
+ CDTP=Z
+ IF(NJ.LE.0) RETURN
+ JL=JF+NJ-1
+ IR=(K-1)*ISTEP
+ DO 3 J=JF,JL
+ 3 CDTP=CDTP+AV((J-1)*ISTEP+I)*AV(J+IR)
+ RETURN
+ END
+CCC
\ No newline at end of file