      PROGRAM CLU
CCC   111130
CCC   SUPPORTS EXPERIMENTAL DIELECTRIC FUNCTIONS ONLY
CCC   NSPH=6; LI=8, LE=8; NXI=200
CCC
CCC   FOR IDFC>=0 ONLY,
CCC   WHEN JXI=JWTM WRITES THE JWTM-TH TRANSITION MATRIX
CCC
      IMPLICIT REAL*8(A-H,O-Z)
CCC   COMMON/C1/VH(NCOU*LITPO),VJ0(NSPH*LMTPO),
CCC  1VYHJ(NCOU*LITPOS),VYJ0(NSPH*LMTPOS),
CCC  2RMI(LI,NSPH),REI(LI,NSPH),W(NLEMT,4),
CCC  3AM0M(NLEMT,NLEMT),FSAS(NSPH),SAS(NSPH,2,2),VINTS(NSPH,16),
CCC  4V3J0(NV3J),SSCS(NSPH),SEXS(NSPH),SABS(NSPH),
CCC  5SQSCS(NSPH),SQEXS(NSPH),SQABS(NSPH),GCSV(NSPH),
CCC  6RXX(NSPH),RYY(NSPH),RZZ(NSPH),ROS(NSPH),RC(NSPH,NSHL),
CCC  7IND3J(LMPO,LM),IOG(NSPH),NSHL(NSPH)
      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
     3FSAC(2,2),SAC(2,2),FSACM(2,2),
     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
     8IND3J(9,8),IOG(6),NSHL(6)
      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
CCC   COMMON/C2/RIS(NHSPO),DLRI(NHSPO),DC0(NSHL-NTL+1),
CCC  1VKT(NSPH),VSZ(NSPH)
CCC   NHSPO=MAX0(NPNT,NPNTTS)*2-1
      COMMON/C2/RIS(999),DLRI(999),DC0(5),VKT(6),VSZ(6)
      COMPLEX*16 RIS,DLRI,DC0,VKT
      COMMON/C3/TFSAS,TSAS(2,2),GCS,SCS,ECS,ACS
      COMPLEX*16 TFSAS,TSAS
      COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
CCC   COMMON/C6/RAC3J(LMTPO)
      COMMON/C6/RAC3J(17)
CCC   COMMON/C9/GIS(NDI,NLEM),GLS(NDI,NLEM),SAM(NDIT,NLEMT)
      COMMON/C9/GIS(480,80),GLS(480,80),SAM(960,160)
      COMPLEX*16 GIS,GLS,SAM
CCC   DIMENSION AM(MXNDM,MXNDM)
CCC   NDIT.LE.MXNDM
      DIMENSION AM(960,960)
      COMPLEX*16 AM
CCC   DIMENSION RCF(NSPH,NSHL)
      DIMENSION RCF(6,8),CMULLR(4,4),CMUL(4,4),CEXTLR(4,4),CEXT(4,4),
     1U(3),US(3),UN(3),UNS(3),UP(3),UPS(3),
     2UNMP(3),UNSMP(3),UPMP(3),UPSMP(3),DUK(3),ARGI(1),ARGS(1)
CCC   DIMENSION ZPV(LM,3,2,2),GAPS(NSPH)
      DIMENSION ZPV(8,3,2,2),GAPS(6),GAPV(3),
     1GAP(3,2),GAPP(3,2),GAPM(3,2),GAPPM(3,2)
      COMPLEX*16 GAPP,GAPPM
CCC   DIMENSION TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH)
      DIMENSION TQSE(2,6),TQSPE(2,6),TQSS(2,6),TQSPS(2,6),
     1TQEV(3),TQSV(3),TQCE(2,3),TQCPE(2,3),TQCS(2,3),TQCPS(2,3)
      COMPLEX*16 TQSPE,TQSPS,TQCPE,TQCPS
CCC   DIMENSION DC0M(NSPH,NSHL-NTL),XIV(NXI)
      DIMENSION DC0M(6,4),XIV(200)
      COMPLEX*16 DC0M
      COMPLEX*16 ARG,S0,S0M,CCSAM
 8010 FORMAT(1H ,16I5)
 8015 FORMAT(1H ,
     1'READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM')
 8050 FORMAT(3(1PD17.8))
 8060 FORMAT(1H ,'READ(IR,*)RXX(I),RYY(I),RZZ(I)')
 8068 FORMAT(1H ,'READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST')
 8069 FORMAT(1H ,'READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST')
 8170 FORMAT(1H ,7(1PD10.3))
 8171 FORMAT(1H ,'READ(IR,*)JWTM')
 6000 FORMAT(1H )
 6050 FORMAT(1H ,' SIZE=',1PD15.7,', REFRACTIVE INDEX=',2(1PD15.7))
 6055 FORMAT(1H ,' SIZE=',1PD15.7)
 6060 FORMAT(1H ,' CIRC ',I2)
 6061 FORMAT(1H ,' CIRC')
 6065 FORMAT(1H ,'  LIN ',I2)
 6066 FORMAT(1H ,'  LIN')
 6068 FORMAT(1H ,4X,'SPHERES; LMX=LI')
 6069 FORMAT(1H ,4X,'SPHERES; LMX=MIN0(LI,LE)')
 6070 FORMAT(1H ,4X,'SPHERE ',I2)
 6080 FORMAT(1H ,4X,'CLUSTER')
 6100 FORMAT(1H ,'----- SCS ----- ABS ----- EXS ----- ALBEDS --')
 6150 FORMAT(1H ,1PD14.7,5(1PD15.7))
 6165 FORMAT(1H ,' UN=(',1PD12.5,',',1PD12.5,',',1PD12.5,')')
 6166 FORMAT(1H ,' UNI=(',1PD12.5,',',1PD12.5,',',1PD12.5,')')
 6167 FORMAT(1H ,' UNS=(',1PD12.5,',',1PD12.5,',',1PD12.5,')')
 6200 FORMAT(1H ,'---- SCS/GS -- ABS/GS -- EXS/GS ---')
 6240 FORMAT(1H ,1X,'FSAS=',2(1PD15.7))
 6242 FORMAT(1H ,1X,'SAS(1,1)=',2(1PD15.7),', SAS(2,1)=',2(1PD15.7))
 6243 FORMAT(1H ,1X,'SAS(1,2)=',2(1PD15.7),', SAS(2,2)=',2(1PD15.7))
 6250 FORMAT(1H ,1X,'FSAT=',2(1PD15.7))
 6252 FORMAT(1H ,1X,'SAT(1,1)=',2(1PD15.7),', SAT(2,1)=',2(1PD15.7))
 6253 FORMAT(1H ,1X,'SAT(1,2)=',2(1PD15.7),', SAT(2,2)=',2(1PD15.7))
 6300 FORMAT(1H ,' STOP IN LUCIN: MATRIX SINGULAR')
 6331 FORMAT(1H ,1X,'MULS') 
 6332 FORMAT(1H ,1X,'MULC') 
 6336 FORMAT(1H ,1X,'MULSLR') 
 6337 FORMAT(1H ,1X,'MULCLR') 
 6350 FORMAT((1H ,7X,4(1PD15.7))) 
 6460 FORMAT(1H ,4X,'CLUSTER (ENSEMBLE AVERAGE, MODE',I2,')')
 6500 FORMAT(1H ,'----- SCC ----- ABC ----- EXC ----- ALBEDC --')
 6550 FORMAT
     1(1H ,' AT ',I1,' LCALC=',I3,' TOO SMALL WITH ARG=',2(1PD15.7))
 6600 FORMAT(1H ,'--- SCC/TGS - ABC/TGS - EXC/TGS ---')
 6650 FORMAT(1H ,' STOP IN HJV')
 6700 FORMAT(1H ,'---- SCCRT --- ABCRT --- EXCRT ----')
 6710 FORMAT(1H ,' REFR. INDEX OF EXTERNAL MEDIUM=',1PD15.7)
 6714 FORMAT(1H ,' VK=',1PD15.7,', XI IS SCALE FACTOR FOR LENGTHS')
 6715 FORMAT(1H ,' VK=',1PD15.7,', XI=',1PD15.7)
 6716 FORMAT(1H ,' XI=',1PD15.7)
 6720 FORMAT(1H ,' TIDG=',1PD10.3,', PIDG=',1PD10.3,
     1', TSDG=',1PD10.3,', PSDG=',1PD10.3)
 6721 FORMAT(1H ,' CFMP=',1PD15.7,', SFMP=',1PD15.7)
 6722 FORMAT(1H ,' CFSP=',1PD15.7,', SFSP=',1PD15.7)
 6740 FORMAT(1H ,' SCAND=',1PD10.3)
 6750 FORMAT(1H ,' STOP IN DME')
 6800 FORMAT(1H ,1X,2('FSAC(',I1,',',I1,')=',2(1PD15.7),3X))
 6810 FORMAT(1H ,1X,2(' SAC(',I1,',',I1,')=',2(1PD15.7),3X))
 6850 FORMAT(1H ,1X,'RE(FSAC(',I1,',',I1,'))/RE(TFSAS)=',1PD15.7,
     1', IM(FSAC(',I1,',',I1,'))/IM(TFSAS)=',1PD15.7)
 6930 FORMAT(1H ,1X,'QSCHU=',1PD15.7,', PSCHU=',1PD15.7,
     1', S0MAG=',1PD15.7)
 6940 FORMAT(1H ,1X,'(RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=',
     11PD15.7)
 6945 FORMAT(1H ,1X,'(IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=',
     11PD15.7)
 6970 FORMAT(1H ,' COSAV=',1PD15.7,', RAPRS=',1PD15.7)
 6971 FORMAT(1H ,' Fx=',1PD15.7,', Fy=',1PD15.7,', Fz=',1PD15.7)
 6973 FORMAT(1H ,' Fl=',1PD15.7,', Fr=',1PD15.7,', Fk=',1PD15.7)
 6974 FORMAT(1H ,' Fk=',1PD15.7)
 6977 FORMAT(1H ,' IPO=',I2,', TQEk=',1PD15.7,', TQSk=',1PD15.7)
 6981 FORMAT
     1(1H ,'  TQEx=',1PD15.7,',  TQEy=',1PD15.7,',  TQEz=',1PD15.7)
 6982 FORMAT
     1(1H ,'  TQSx=',1PD15.7,',  TQSy=',1PD15.7,',  TQSz=',1PD15.7)
 6983 FORMAT
     1(1H ,'  TQEl=',1PD15.7,',  TQEr=',1PD15.7,',  TQEk=',1PD15.7)
 6984 FORMAT
     1(1H ,'  TQSl=',1PD15.7,',  TQSr=',1PD15.7,',  TQSk=',1PD15.7)
 6991 FORMAT('========== JXI =',I3,' ====================')
 6992 FORMAT('********** JTH =',I3,', JPH =',I3,
     1', JTHS =',I3,', JPHS =',I3,' ********************')
      IR=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   MXNDM IS FIRST NUMERICAL DECLARED DIMENSION OF AM;
CCC   NDIT=2*NDI; NDI=NSPH*NLIM
CCC
CCC
CCC   VK*ROS(I) IS SIZE PARAMETER in vacuo
CCC
CCC
CCC   OUTPUTS (NORMALIZED SCATTERING AMPLITUDE)*(2*PIG/(VK*VK)) IN m**3;
CCC   OUTPUTS FORCE*(c/I0) IN m**2; 
CCC   OUTPUTS (EXT_TO_TORQUE) AND (SCA_TO_TORQUE) CONTRIBUTIONS 
CCC   *(c/I0)*(VK*EXRI)**3; (TORQUE OUTPUTS)/(EXDC*VK*VK) ARE IN m**2
CCC
CCC
CCC   ON A LINEAR BASIS,
CCC   POLARIZATION (1), (1,ni), (ns,1) OR -1 IS PARALLE(l),
CCC   POLARIZATION (2), (2,ni), (ns,2) OR  1 IS PERPENDICULA(r);
CCC   ON A CIRCULAR BASIS,
CCC   POLARIZATION (1), (1,ni), (ns,1) OR -1 IS (R)ight,
CCC   POLARIZATION (2), (2,ni), (ns,2) OR  1 IS (L)eft;
CCC   ns=SCATWAVE POLARIZATION, ni=INCWAVE POLARIZATION
CCC
CCC
CCC   NV3J=(LM*LMPO*(LMT+7))/6;
CCC   LM=MAX0(LI,LE)
CCC
CCC
CCC   IAVM=0 --> RANDOM AVERAGES FOR AMPLITUDE AND CROSS SECTIONS;
CCC   IAVM=1 --> SAME AS IAVM=0, PLUS RANDOM AVERAGES FOR INTENSITY
CCC
CCC
CCC   ISAM<0 --> ARBITRARILY VARYING SCATTERING PLANE,
CCC              AS IT IS UNDEFINED FOR FORWARD OR BACKWARD SCATTERING,
CCC              MERIDIONAL PLANE FOR INCIDENCE ANGLES IS ASSUMED;
CCC   ISAM=0 --> MERIDIONAL PLANES
CCC              (I. E. REFERENCE PLANES ACROSS Z AXIS)
CCC              AS DETERMINED BY INCIDENCE AND OBSERVATION ANGLES;
CCC   ISAM=1 --> SCATTERING PLANE ACROSS THE Z AXIS; SCATTERING
CCC              AS A FUNCTION OF INCIDENCE AND OF OBSERVATION ANGLES;
CCC              NEEDS ONLY TH, PH, THS; THSLST GT 180 IF NECESSARY;
CCC   ISAM>1 --> SCATTERING PLANE ACROSS THE Z AXIS; SCATTERING
CCC              AS A FUNCTION OF INCIDENCE ANGLES
CCC              FOR A FIXED SCATTERING ANGLE THSCA=THS-TH;
CCC              NEEDS ONLY TH, PH, THS
CCC
      OPEN(IR,FILE='../../test_data/cluster/DCLU',STATUS='OLD')
      READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM
      DO 101 I=1,NSPH
  101 READ(IR,*)RXX(I),RYY(I),RZZ(I)
      READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST
      READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST
      READ(IR,*)JWTM
      CLOSE(IR)
      OPEN(IW,FILE='OCLU',STATUS='UNKNOWN')
      WRITE(IW,8015)
      WRITE(IW,8010)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM
      WRITE(IW,8060)
      WRITE(IW,8050)(RXX(I),RYY(I),RZZ(I),I=1,NSPH)
      WRITE(IW,8068)
      WRITE(IW,8170)TH,THSTP,THLST,THS,THSSTP,THSLST
      WRITE(IW,8069)
      WRITE(IW,8170)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST
      WRITE(IW,8171)
      WRITE(IW,8010)JWTM
      OPEN(ITIN,FILE='TEDF',FORM='UNFORMATTED',STATUS='OLD')
      READ(ITIN)NSPHT
      IF(NSPHT.NE.NSPH)GO TO 497
      READ(ITIN)(IOG(I),I=1,NSPH)
      READ(ITIN)EXDC,WP,XIP,IDFC,NXI
      READ(ITIN)(XIV(I),I=1,NXI)
      DO 103 I=1,NSPH
      IF(IOG(I).LT.I)GO TO 103
      READ(ITIN)NSHL(I),ROS(I)
      NSH=NSHL(I)
      READ(ITIN)(RCF(I,NS),NS=1,NSH)
  103 CONTINUE
      WRITE(IW,*)' READ(ITIN)NSPHT'
      WRITE(IW,*)' READ(ITIN)(IOG(I),I=1,NSPH)'
      WRITE(IW,*)' READ(ITIN)EXDC,WP,XIP,IDFC,NXI'
      WRITE(IW,*)' READ(ITIN)(XIV(I),I=1,NXI)'
      WRITE(IW,*)' READ(ITIN)NSHL(I),ROS(I)'
      WRITE(IW,*)' READ(ITIN)(RCF(I,NS),NS=1,NSH)'
      WRITE(IW,6000)
      C0=0.0D0
      SML=1.0D-3
      IF(THSTP.EQ.C0)NTH=0
      IF(THSTP.NE.C0)NTH=IDINT((THLST-TH)/THSTP+SML)
      NTH=NTH+1
      IF(PHSTP.EQ.C0)NPH=0
      IF(PHSTP.NE.C0)NPH=IDINT((PHLST-PH)/PHSTP+SML)
      NPH=NPH+1
      IF(ISAM.LE.1)GO TO 110 
      NTHS=1
      THSCA=THS-TH
      GO TO 111
  110 IF(THSSTP.EQ.C0)NTHS=0
      IF(THSSTP.NE.C0)NTHS=IDINT((THSLST-THS)/THSSTP+SML)
      NTHS=NTHS+1
      IF(ISAM.LT.1)GO TO 112
  111 NPHS=1
      GO TO 113
  112 IF(PHSSTP.EQ.C0)NPHS=0
      IF(PHSSTP.NE.C0)NPHS=IDINT((PHSLST-PHS)/PHSSTP+SML)
      NPHS=NPHS+1
  113 NK=NTH*NPH
      NKS=NTHS*NPHS
      NKKS=NK*NKS
      TH1=TH
      PH1=PH
      THS1=THS
      PHS1=PHS
      CALL STR(RCF)
      LM=MAX0(LI,LE)
      CALL THDPS(LM,ZPV)
      EXRI=DSQRT(EXDC)
      WRITE(IW,6710)EXRI
      OPEN(IT,FILE='TPPOAN',FORM='UNFORMATTED',STATUS='UNKNOWN')
      WRITE(IT)IAVM,ISAM,INPOL
      WRITE(IT)NXI,NTH,NPH,NTHS,NPHS

C
      NLEMT=NLEM+NLEM
C

      WN=WP/3.0D08
      SQSFI=1.0D0
      IF(IDFC.GE.0)GO TO 118
      VK=XIP*WN
      WRITE(IW,6714)VK
      WRITE(IW,6000)
  118 DO 488 JXI=1,NXI
      JAW=1
      WRITE(IW,6991)JXI
      XI=XIV(JXI)
      IF(IDFC.LT.0)GO TO 119
      VK=XI*WN
      VKARG=VK
      WRITE(IW,6715)VK,XI
      GO TO 120
  119 VKARG=XI*VK
      SQSFI=1.0D0/(XI*XI)
      WRITE(IW,6716)XI
 120  CALL HJV(EXRI,VKARG,JER,LCALC,ARG)
C     DEBUG CODE
C     IF(JXI.NE.1)GOTO 199
C     PRINT *,"DEBUG: witing VJ0"
C     CALL LOGMAT(VJ0)
C 199 PRINT *,"DEBUG: already wrote VJ0"
C     END DEBUG
      IF(JER.EQ.0)GO TO 122
      WRITE(IW,6650)
      GO TO 490
  122 DO 132 I=1,NSPH
      IOGI=IOG(I)
      IF(IOGI.EQ.I)GO TO 125
      DO 123 L=1,LI
      RMI(L,I)=RMI(L,IOGI)
  123 REI(L,I)=REI(L,IOGI)
      GO TO 132
  125 NSH=NSHL(I)
      ICI=(NSH+1)/2
      IF(IDFC.NE.0)GO TO 127
      READ(ITIN)(DC0(IC),IC=1,ICI)
      GO TO 130
  127 IF(JXI.EQ.1)READ(ITIN)(DC0M(I,IC),IC=1,ICI)
      DO 128 IC=1,ICI
  128 DC0(IC)=DC0M(I,IC)
  130 IF(MOD(NSH,2).EQ.0)DC0(ICI+1)=EXDC
      CALL DME(LI,I,NPNT,NPNTTS,VKARG,EXDC,EXRI,JER,LCALC,ARG)
C      DO 199 IDL=1,8
C      PRINT *,"DEBUG: RMI(",IDL,",",I,") =",RMI(IDL,I)
C     199  PRINT *,"DEBUG: REI(",IDL,",",I,") =",REI(IDL,I)
C      PRINT *,"DEBUG: DME - jer =",JER,", lcalc =",LCALC,
C     1", arg =",ARG
      IF(JER.EQ.0)GO TO 132
      WRITE(IW,6750)
      GO TO 490
 132  CONTINUE
      CALL CMS(AM)
      NDIT=2*NSPH*NLIM
      CALL LUCIN(AM,MXNDM,NDIT,JER)
      IF(JER.EQ.1)GO TO 495
      CALL ZTM(AM)

C
      IF(IDFC.LT.0)GO TO 156
      IF(JXI.NE.JWTM)GO TO 156
      IS=1
      ITS=8
      OPEN(ITS,FILE='TTMS',FORM='UNFORMATTED',STATUS='UNKNOWN')
      WRITE(ITS)IS,LM
      WRITE(ITS)VK,EXRI
      WRITE(ITS)((AM0M(I,J),J=1,NLEMT),I=1,NLEMT)
      CLOSE(ITS)
  156 CONTINUE
C

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