      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=25
      IW=26
      IT=27
      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='../../test_data/sphere/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
      PRINT *,"INFO: done scale iteration ",JXI
  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)
      PRINT *,"INFO: output written to OSPH."
      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
      IF(ZP.EQ.C0) GOTO 15
C      PRINT 9010,L-1,ZP
C      PRINT 9011,L-1,ZP
 9010 FORMAT("DEBUG: zpv[",I1,"][1][0][1] = ",E15.3)
 9011 FORMAT("DEBUG: zpv[",I1,"][1][1][0] = ",E15.3)
   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
      IF(ZP.EQ.C0) GOTO 20
C      PRINT 9012,L-1,ZP
C      PRINT 9013,L-1,ZP
 9012 FORMAT("DEBUG: zpv[",I1,"][0][0][0] = ",E15.3)
 9013 FORMAT("DEBUG: zpv[",I1,"][0][1][1] = ",E15.3)
   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
      IF(ZP.EQ.C0) GOTO 25
C      PRINT 9014,L-1,ZP
C      PRINT 9015,L-1,ZP
 9014 FORMAT("DEBUG: zpv[",I1,"][2][0][0] = ",E15.3)
 9015 FORMAT("DEBUG: zpv[",I1,"][2][1][1] = ",E15.3)
   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
