Skip to content
Snippets Groups Projects
Commit dfd4aaa6 authored by Giacomo Mulas's avatar Giacomo Mulas
Browse files

Upload New File

parent bff81c0b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 1493 additions and 0 deletions
sphere/sph.f 0 → 100644
+ 1493
− 0
View file @ dfd4aaa6
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment