From dfd4aaa603774c72eea569cbc6462f318c952e17 Mon Sep 17 00:00:00 2001
From: Giacomo Mulas <giacomo.mulas@inaf.it>
Date: Mon, 4 Sep 2023 12:19:32 +0000
Subject: [PATCH] Upload New File
---
sphere/sph.f | 1493 ++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 1493 insertions(+)
create mode 100644 sphere/sph.f
diff --git a/sphere/sph.f b/sphere/sph.f
new file mode 100644
index 00000000..f7aafd52
--- /dev/null
+++ b/sphere/sph.f
@@ -0,0 +1,1493 @@
+ PROGRAM SPH
+CCC 140531
+CCC SUPPORTS EXPERIMENTAL DIELECTRIC FUNCTIONS ONLY
+CCC NSPH=2; LM=40; NXI=200
+CCC
+CCC FOR IDFC>=0 AND NSPH=1 ONLY,
+CCC WHEN JXI=JWTM WRITES THE JWTM-TH TRANSITION MATRIX
+CCC
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC COMMON/C1/RMI(LM,NSPH),REI(LM,NSPH),W(NLMMT,4),
+CCC 1FSAS(NSPH),SAS(NSPH,2,2),VINTS(NSPH,16),
+CCC 2SSCS(NSPH),SEXS(NSPH),SABS(NSPH),
+CCC 3SQSCS(NSPH),SQEXS(NSPH),SQABS(NSPH),GCSV(NSPH),
+CCC 4ROS(NSPH),RC(NSPH,NSHL),IOG(NSPH),NSHL(NSPH)
+ COMMON/C1/RMI(40,2),REI(40,2),W(3360,4),
+ 1FSAS(2),SAS(2,2,2),VINTS(2,16),
+ 2SSCS(2),SEXS(2),SABS(2),SQSCS(2),SQEXS(2),SQABS(2),
+ 3GCSV(2),RXX(1),RYY(1),RZZ(1),ROS(2),RC(2,8),
+ 4IOG(2),NSHL(2)
+ COMPLEX*16 RMI,REI,W,FSAS,SAS,VINTS
+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(2),VSZ(2)
+ COMPLEX*16 RIS,DLRI,DC0,VKT
+CCC DIMENSION RCF(NSPH,NSHL)
+ DIMENSION RCF(2,8),CMULLR(4,4),CMUL(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)
+ DIMENSION VINT(16)
+ COMPLEX*16 VINT
+CCC DIMENSION ZPV(LM,3,2,2),GAPS(NSPH)
+ DIMENSION ZPV(40,3,2,2),GAPS(2)
+CCC DIMENSION TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH)
+ DIMENSION TQSE(2,2),TQSPE(2,2),TQSS(2,2),TQSPS(2,2)
+ COMPLEX*16 TQSPE,TQSPS
+CCC DIMENSION DC0M(NSPH,NSHL-NTL),XIV(NXI)
+ DIMENSION DC0M(2,4),XIV(200)
+ COMPLEX*16 DC0M
+ COMPLEX*16 ARG,S0,TFSAS
+ 8010 FORMAT(1H ,16I5)
+ 8015 FORMAT(1H ,'READ(IR,*)NSPH,LM,INPOL,NPNT,NPNTTS,ISAM')
+ 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)
+ 6061 FORMAT(1H ,' CIRC')
+ 6066 FORMAT(1H ,' LIN')
+ 6070 FORMAT(1H ,4X,'SPHERE ',I2)
+ 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))
+ 6331 FORMAT(1H ,1X,'MULS')
+ 6336 FORMAT(1H ,1X,'MULSLR')
+ 6350 FORMAT((1H ,7X,4(1PD15.7)))
+ 6550 FORMAT
+ 1(1H ,' AT ',I1,' LCALC=',I3,' TOO SMALL WITH ARG=',2(1PD15.7))
+ 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')
+ 6930 FORMAT(1H ,1X,'QSCHU=',1PD15.7,', PSCHU=',1PD15.7,
+ 1', S0MAG=',1PD15.7)
+ 6970 FORMAT(1H ,' COSAV=',1PD15.7,', RAPRS=',1PD15.7)
+ 6971 FORMAT(1H ,' Fx=',1PD15.7,', Fy=',1PD15.7,', Fz=',1PD15.7)
+ 6977 FORMAT(1H ,' IPO=',I2,', TQEk=',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 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 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='DSPH',STATUS='OLD')
+ READ(IR,*)NSPH,LM,INPOL,NPNT,NPNTTS,ISAM
+ READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST
+ READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST
+ READ(IR,*)JWTM
+ CLOSE(IR)
+ OPEN(IW,FILE='OSPH',STATUS='UNKNOWN')
+ WRITE(IW,8015)
+ WRITE(IW,8010)NSPH,LM,INPOL,NPNT,NPNTTS,ISAM
+ 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
+ PIGH=DACOS(C0)
+ PIG=PIGH*2.0D0
+ GCS=C0
+ DO 116 I=1,NSPH
+ IOGI=IOG(I)
+ IF(IOGI.LT.I)GO TO 116
+ GCSS=PIG*ROS(I)*ROS(I)
+ GCSV(I)=GCSS
+ NSH=NSHL(I)
+ DO 115 J=1,NSH
+ 115 RC(I,J)=RCF(I,J)*ROS(I)
+ 116 GCS=GCS+GCSV(IOGI)
+ CALL THDPS(LM,ZPV)
+ EXRI=DSQRT(EXDC)
+ WRITE(IW,6710)EXRI
+ OPEN(IT,FILE='TPPOAN',FORM='UNFORMATTED',STATUS='UNKNOWN')
+ IMODE=10
+ WRITE(IT)IMODE,ISAM,INPOL
+ WRITE(IT)NXI,NTH,NPH,NTHS,NPHS
+ WRITE(IT)NSPH
+ WRITE(IT)(IOG(NS),NS=1,NSPH)
+ IF(INPOL.EQ.0)WRITE(IW,6066)
+ IF(INPOL.NE.0)WRITE(IW,6061)
+ WRITE(IW,6000)
+ 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
+ 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 WRITE(IT)VK
+ DO 132 I=1,NSPH
+ IOGI=IOG(I)
+ IF(IOGI.EQ.I)GO TO 125
+ DO 123 L=1,LM
+ 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(LM,I,NPNT,NPNTTS,VKARG,EXDC,EXRI,JER,LCALC,ARG)
+ IF(JER.EQ.0)GO TO 132
+ WRITE(IW,6750)
+ GO TO 490
+ 132 CONTINUE
+
+C
+ IF(IDFC.LT.0)GO TO 156
+ IF(NSPH.GT.1)GO TO 156
+ IF(JXI.NE.JWTM)GO TO 156
+ IS=1111
+ ITS=8
+ OPEN(ITS,FILE='TTMS',FORM='UNFORMATTED',STATUS='UNKNOWN')
+ WRITE(ITS)IS,LM
+ WRITE(ITS)VK,EXRI
+ WRITE(ITS)(-1.0D0/RMI(I,1),-1.0D0/REI(I,1),I=1,LM)
+ WRITE(ITS)ROS(1)
+ CLOSE(ITS)
+ 156 CONTINUE
+C
+
+ CS0=0.25D0*VK*VK*VK/PIGH
+ CALL SSCR0(TFSAS,NSPH,LM,VK,EXRI)
+ 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)
+ 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)
+ WRITE(IT)TQSE(1,I),TQSS(1,I),TQSPE(1,I),TQSPS(1,I)
+ WRITE(IT)TQSE(2,I),TQSS(2,I),TQSPE(2,I),TQSPS(2,I)
+ 170 CONTINUE
+ IF(NSPH.EQ.1)GO TO 172
+ 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
+ 172 TH=TH1
+ DO 486 JTH=1,NTH
+ PH=PH1
+ DO 484 JPH=1,NPH
+ IF((NK.EQ.1).AND.(JXI.GT.1))GO TO 182
+ CALL UPVMP(TH,PH,0,COST,SINT,COSP,SINP,U,UPMP,UNMP)
+ 182 IF(ISAM.LT.0)GO TO 184
+ CALL WMAMP
+ 1(0,COST,SINT,COSP,SINP,INPOL,LM,0,NSPH,ARGI,U,UPMP,UNMP)
+ DO 183 I=1,NSPH
+ RAPR=SEXS(I)-GAPS(I)
+ FRX=RAPR*U(1)
+ FRY=RAPR*U(2)
+ FRZ=RAPR*U(3)
+ 183 CONTINUE
+ 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,LM,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,LM,0,NSPH,ARGI,ARGS)
+ GO TO 194
+ 192 DO 193 I=1,3
+ UN(I)=UNMP(I)
+ 193 UNS(I)=UNSMP(I)
+ 194 IF(ISAM.LT.0)JW=1
+ WRITE(IT)TH,PH,THS,PHS
+ WRITE(IT)SCAN
+ IF(JW.EQ.0)GO TO 212
+ JW=0
+ WRITE(IT)U(1),U(2),U(3)
+ 212 CONTINUE
+ 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 224
+ 214 WRITE(IW,6165)UN(1),UN(2),UN(3)
+ 224 CALL SSCR2(NSPH,LM,VK,EXRI)
+ DO 226 NS=1,NSPH
+ IF(IOG(NS).LT.NS)GO TO 226
+ WRITE(IW,6070)NS
+ WRITE(IW,6242)SAS(NS,1,1),SAS(NS,2,1)
+ WRITE(IW,6243)SAS(NS,1,2),SAS(NS,2,2)
+ IF((JTHS.EQ.1).AND.(JPHS.EQ.1))WRITE(IW,6971)FRX,FRY,FRZ
+ DO 225 I=1,16
+ 225 VINT(I)=VINTS(NS,I)
+ CALL MMULC(VINT,CMULLR,CMUL)
+ WRITE(IW,6331)
+ WRITE(IW,6350)((CMUL(I,J),J=1,4),I=1,4)
+ WRITE(IW,6336)
+ WRITE(IW,6350)((CMULLR(I,J),J=1,4),I=1,4)
+ WRITE(IT)(VINT(I),I=1,16)
+ WRITE(IT)((CMUL(I,J),J=1,4),I=1,4)
+ 226 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
+ 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(2),RMI(40,2),REI(40,2),
+ 1TQSE(2,2),TQSPE(2,2),TQSS(2,2),TQSPS(2,2)
+ 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 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 SSCR0(TFSAS,NSPH,LM,VK,EXRI)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C1/RMI(40,2),REI(40,2),W(3360,4),
+ 1FSAS(2),SAS(2,2,2),VINTS(2,16),
+ 2SSCS(2),SEXS(2),SABS(2),SQSCS(2),SQEXS(2),SQABS(2),
+ 3GCSV(2),RXX(1),RYY(1),RZZ(1),ROS(2),RC(2,8),
+ 4IOG(2),NSHL(2)
+ COMPLEX*16 RMI,REI,W,FSAS,SAS,VINTS
+ COMPLEX*16 TFSAS,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)
+ TFSAS=CC0
+ DO 12 I=1,NSPH
+ IOGI=IOG(I)
+ IF(IOGI.LT.I)GO TO 12
+ SUMS=C0
+ SUM21=CC0
+ DO 10 L=1,LM
+ 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 TFSAS=TFSAS+FSAS(IOGI)
+ RETURN
+ END
+ SUBROUTINE SSCR2(NSPH,LM,VK,EXRI)
+ IMPLICIT REAL*8(A-H,O-Z)
+ COMMON/C1/RMI(40,2),REI(40,2),W(3360,4),
+ 1FSAS(2),SAS(2,2,2),VINTS(2,16),
+ 2SSCS(2),SEXS(2),SABS(2),SQSCS(2),SQEXS(2),SQABS(2),
+ 3GCSV(2),RXX(1),RYY(1),RZZ(1),ROS(2),RC(2,8),
+ 4IOG(2),NSHL(2)
+ COMPLEX*16 RMI,REI,W,FSAS,SAS,VINTS
+ COMPLEX*16 S11,S21,S12,S22,RM,RE,CC0,CSAM,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)
+ NLMM=LM*(LM+2)
+ DO 14 I=1,NSPH
+ IOGI=IOG(I)
+ IF(IOGI.LT.I)GO TO 14
+ K=0
+ S11=CC0
+ S21=CC0
+ S12=CC0
+ S22=CC0
+ DO 10 L=1,LM
+ 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+NLMM
+ 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
+ 14 CONTINUE
+ 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
+ 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(40,3,2,2),IOG(2),RMI(40,2),REI(40,2),GAPS(2)
+ 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 THDPS(LM,ZPV)
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC DIMENSION ZPV(LM,3,2,2)
+ DIMENSION ZPV(40,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(1682),U(3),UP(3),UN(3),ARG(1)
+ COMPLEX*16 YLM
+ COMMON/C1/RMI(40,2),REI(40,2),W(3360,4),
+ 1FSAS(2),SAS(2,2,2),VINTS(2,16),
+ 2SSCS(2),SEXS(2),SABS(2),SQSCS(2),SQEXS(2),SQABS(2),
+ 3GCSV(2),RXX(1),RYY(1),RZZ(1),ROS(2),RC(2,8),
+ 4IOG(2),NSHL(2)
+ COMPLEX*16 RMI,REI,W,FSAS,SAS,VINTS
+ 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(1682),ARGI(1),ARGS(1),
+ 1U(3),UP(3),UN(3),US(3),UPS(3),UNS(3)
+ COMPLEX*16 YLM
+ COMMON/C1/RMI(40,2),REI(40,2),W(3360,4),
+ 1FSAS(2),SAS(2,2,2),VINTS(2,16),
+ 2SSCS(2),SEXS(2),SABS(2),SQSCS(2),SQEXS(2),SQABS(2),
+ 3GCSV(2),RXX(1),RYY(1),RZZ(1),ROS(2),RC(2,8),
+ 4IOG(2),NSHL(2)
+ COMPLEX*16 RMI,REI,W,FSAS,SAS,VINTS
+ 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 PWMA(UP,UN,YLM,INPOL,LW,ISQ)
+ IMPLICIT REAL*8(A-H,O-Z)
+CCC DIMENSION YLM(NLWM+2)
+ DIMENSION YLM(1682),UP(3),UN(3)
+ COMPLEX*16 YLM
+ COMMON/C1/RMI(40,2),REI(40,2),W(3360,4),
+ 1FSAS(2),SAS(2,2,2),VINTS(2,16),
+ 2SSCS(2),SEXS(2),SABS(2),SQSCS(2),SQEXS(2),SQABS(2),
+ 3GCSV(2),RXX(1),RYY(1),RZZ(1),ROS(2),RC(2,8),
+ 4IOG(2),NSHL(2)
+ COMPLEX*16 RMI,REI,W,FSAS,SAS,VINTS
+ 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 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(42),RFJ(42),RFN(42),
+ 1FBI(42),FB(42),FN(42),RMF(40),DRMF(40),REF(40),DREF(40)
+ COMPLEX*16 CFJ,FBI,FB,FN,RMF,DRMF,REF,DREF
+ COMMON/C1/RMI(40,2),REI(40,2),W(3360,4),
+ 1FSAS(2),SAS(2,2,2),VINTS(2,16),
+ 2SSCS(2),SEXS(2),SABS(2),SQSCS(2),SQEXS(2),SQABS(2),
+ 3GCSV(2),RXX(1),RYY(1),RZZ(1),ROS(2),RC(2,8),
+ 4IOG(2),NSHL(2)
+ COMPLEX*16 RMI,REI,W,FSAS,SAS,VINTS
+ COMMON/C2/RIS(999),DLRI(999),DC0(5),VKT(2),VSZ(2)
+ 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(2),VSZ(2)
+ 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/RMI(40,2),REI(40,2),W(3360,4),
+ 1FSAS(2),SAS(2,2,2),VINTS(2,16),
+ 2SSCS(2),SEXS(2),SABS(2),SQSCS(2),SQEXS(2),SQABS(2),
+ 3GCSV(2),RXX(1),RYY(1),RZZ(1),ROS(2),RC(2,8),
+ 4IOG(2),NSHL(2)
+ COMPLEX*16 RMI,REI,W,FSAS,SAS,VINTS
+ COMMON/C2/RIS(999),DLRI(999),DC0(5),VKT(2),VSZ(2)
+ 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 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=LM
+ DIMENSION SINRMP(40),COSRMP(40),PLEGN(861),YLM(1681)
+ 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
+CCC
\ No newline at end of file
--
GitLab