***********************************************************************
      PROGRAM MOMDIS
* Calculates cross sections and momentum distributions in knockout 
* reactions at intemediate energies (30 MeV < E_lab/nucl < 2200 MeV).
* Based on the article:
* C.A. Bertulani and P.G. Hansen, Phys. Rev.  C70, 034609 (2004) 
*
* A `readme' file with a description of how to prepare an input is 
* attached at the end of this code. 
*
* Written by C.A. Bertulani and A. Gade,  2006
*======================================================================
*
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 IV(0:40),FAK(0:40)
      COMPLEX*16 CI 
      COMMON/A1/ BETA,GAMMA
      COMMON /SYST/ APROJ,ZPROJ,ATARG,ZTARG,ACORE,ZCORE,ECA,S2C
      COMMON /COUL/ PCM,PCM2,ETA,ETA2
      COMMON /PMOM/ PMAX
      COMMON /CONST1/ PI,FOUR_PI,HBARC,AMU,FINE,E2
      COMMON /MASS/ AVMASS,ACMASS
      COMMON /RATIOS/ RATIOC,RATIOV
      COMMON /CONST2/ CI
      COMMON /ML/ MLML
      COMMON /GF1 / IV
      COMMON /GF2 / FAK
      COMMON /MLC/ CLM
      COMMON /CUT/ BMAX,RMAX
*
      CHARACTER*12 FILENA
*
      PRINT*,' Enter name of INPUT file'
      READ(*,*)FILENA
*
      OPEN(UNIT=12,FILE=FILENA,STATUS='UNKNOWN')
*
* enter reaction parameters
*--------------------------
      CALL SKIPCOM(12)
      READ(12,*)ZPROJ,APROJ,ZTARG,ATARG,ECA
      CALL SKIPCOM(12)
      READ(12,*)ZCORE,ACORE
      CALL SKIPCOM(12)
      READ(12,*)S2C
*
      PRINT*,' '
      PRINT*,' --------------  Program MOMDIS     ----------------' 
      PRINT*,' '
      PRINT*,'    Calculates momentum distributions for stripping' 
      PRINT*,'    and diffraction dissociation reactions'
      PRINT*,' '
      PRINT*,' ---------------------------------------------------' 
      PRINT*,' '
*
* famous constants 
*-----------------
      CI      = CMPLX(0.,1.)
      PI      = 3.14159265359D0
      FOUR_PI = 4.D0*PI
      AMU     = 931.49386D0
      HBARC   = 197.3270533D0
      FINE    = 1./137.036D0
      E2      = HBARC*FINE
*
* variables used in the Coulomb phase of the core+target scattering
*------------------------------------------------------------------
      APMASS = EXCESS(ZPROJ,APROJ,1.D0) + APROJ*AMU - ZPROJ*0.510999
      ATMASS = EXCESS(ZTARG,ATARG,1.D0) + ATARG*AMU - ZTARG*0.510999
      REDM   = APMASS*ATMASS/(APMASS+ATMASS)/AMU
     	ECM    = ELAB*REDM/APROJ
      ELAB   = ECA*APROJ
      PLAB   = SQRT(2.*APMASS*ELAB)/HBARC
      GAMMA  = (ECA+AMU)/AMU
      BETA   = SQRT(1.-1./GAMMA/GAMMA)
      ETA    = ZCORE*ZTARG*FINE/BETA
      ZVAL   = ZPROJ-ZCORE
      AVAL   = APROJ-ACORE
      ETA2   = ZVAL*ZTARG*FINE/BETA
      IF(ZCORE.EQ.-1.) THEN
        ETA  = 0.
        ETA2 = 0.
      ENDIF
      ACMASS = EXCESS(ZCORE,ACORE,1.D0) + ACORE*AMU - ZCORE*0.510999
      AVMASS = EXCESS(ZVAL,AVAL,1.D0) + AVAL*AMU - ZVAL*0.510999
      ARED   = ACMASS*ATMASS/(ACMASS+ATMASS)
      PCM    = ARED/ACMASS*PLAB
      PCM2   = PCM*(ACMASS+ATMASS)/(AVMASS+ATMASS)
* mass ratios, integration limits and cutoffs
*--------------------------------------------
      RATIOC = AVMASS/APMASS 
      RATIOV = ACMASS/APMASS 
      BMAX   = 20.
      RMAX  =  30.
*
* wanna work?
*------------	      
      PRINT*,' Enter:'
      PRINT*,' 1 for energy and wavefunction of bound states'
      PRINT*,' 2 for eikonal t-rho-rho S-matrices'
      PRINT*,' 3 for momentum distributions'
      READ(*,*) IOP
*
      IF(IOP.EQ.1) CALL EIGEN
      IF(IOP.EQ.2) THEN
* want S-matrix for core+target or nucleon+target?
*-------------------------------------------------
        CALL SKIPCOM(12)
        READ(12,*)ISMAT
        ANN=APROJ-ACORE
        ZNN=ZPROJ-ZCORE
        IF(ISMAT.EQ.0)CALL SMAT(ANN,ZNN,ATARG,ZTARG,ECA)
        IF(ISMAT.EQ.1)CALL SMAT(ACORE,ZCORE,ATARG,ZTARG,ECA)
      ENDIF
*
* for momentum distributions
*---------------------------
      IF(IOP.EQ.3) THEN
        PRINT*,' '
        PRINT*,' Enter:'
        PRINT*,'       1 for d sigma / d p_z'
        PRINT*,'       2 for d^2 sigma / d p_t^2'
        PRINT*,'       3 for d sigma / d p_y'
        PRINT*,'       4 for d^2 sigma / d p_z d p_t'
        PRINT*,'       5 for core + target elastic cross section'
        READ(*,*)IOPT
        PRINT*,' '
        IF(IOPT.NE.5) THEN
          PRINT*,' Enter angular momentum L,M of the wavefunction'
          READ(*,*)L,M
        MLML=2
        IF(M.EQ.0)MLML=1
        ENDIF
*
* call routine for signs and factorials
*--------------------------------------
      CALL GFV(40)
*
* CLM coefficients
*-----------------
        MM=IABS(M)
        CLM=IV(MM)*SQRT((2.*L+1.)/4./PI*FAK(L-MM)/FAK(L+MM))
        PRINT*,' '
*
* maximum momentum (in MeV/c)
*----------------------------
        PMAX=400.
* starting integration parameters
*---------------------------------
        NB=110
        NR=110	       
        NK=110
*
        CALL CROSS(L,M,IOPT,NR,NB,NK)
        PRINT*,' '
        PRINT*,' DONE: Now enter any number to stop program'
        READ(*,*)DUM
      ENDIF
*
      STOP
      END
*
*======================================================================
*
      SUBROUTINE EIGEN
*
* finds eigenvalue and eigenfunction of a particle with charge Z1,
* in a WS + spin-orbit + Coulomb potential with a nucleus with charge Z2.
*
*   input parameters
*   ================
* 
* NPNTS  = number of points in the radial mesh
* RMAX   = maximum radius size
* N_0    = nodes of the wave function 
* L0     = orbital angular momentum of single particle state
* J0     = total angular momentum of single particle state 
* V0     = depth of central potential 
* VS0    = depth of spin-orbit potential 
* R0     = radius of the central potential
* RC     = radius of Coulomb potential 
* RS0    = radius of the spin-orbit potential
* AAS    = diffuseness of the spin-orbit potential 
* Z1, Z2 = charges of the nuclei
* A1, A2 = masses of the nuclei (in amu) 
*
      IMPLICIT REAL*8 (A-H,O-Z)
*
* NLG = maximum number of radial coordinate points 
      PARAMETER(NLG = 50000)
*
      CHARACTER ANAME1*5,ANAME2*5
      CHARACTER*12 FILENA
*
      COMMON /POTENTIAL/ U(NLG),WSO(NLG),UE(NLG)
     &                   ,VC(NLG)
      COMMON /WAVE_FUNCTIONS/ EIGWF(NLG)
      COMMON /PAR/  DR,DR2,PI,FOUR_PI,NPNTS
      COMMON /POT/ V0,VS0,R0,RS0,RC,AA,AAS,Z1Z2,ARED
      COMMON /SYST/ APROJ,ZPROJ,ATARG,ZTARG,ACORE,ZCORE,ECA,S2C
      COMMON /MASS/ AVMASS,ACMASS
*
      OPEN(UNIT=13,FILE='EIGEN.TXT', STATUS='UNKNOWN')
*
      PRINT*,' Enter name of OUTPUT file'
      READ(*,*)FILENA
*
      OPEN(UNIT=16,FILE=FILENA,STATUS='UNKNOWN')
*
      WRITE(6,*)' '
      WRITE(13,*)'** Bound-state in cent. + Coul. + spin-orbit pot **'
      WRITE(13,*)' '
*
* nodes of the wave function, total angular momentum, orbital ang. mom.
*----------------------------------------------------------------------
      CALL SKIPCOM(12)
      READ(12,*) N_0,AJ0,L0
*
* mesh size, no. of points in the radial mesh, max. radius size
*--------------------------------------------------------------
      NPNTS = 50000 
      RMAX  = 250.   ! large R_max for occasional 'halo' states
*
      DR=RMAX/NPNTS
*
* potential parameters  
*---------------------
      CALL SKIPCOM(12)
      READ(12,*)V0,R0,AA,VS0,RS0,AAS,RC
*
* masses and charges
*-------------------
      Z1 = ZPROJ-ZCORE
	A1 = APROJ-ACORE
	Z2 = ZCORE
	A2 = ACORE
*
      CLOSE(UNIT=12)
*
      WRITE(13,*)'Potential parameters'
      WRITE(13,*)' V0            R0          AA'
      WRITE(13,20)V0,R0,AA
      WRITE(13,*)' VS0           RS0         AAS          RC'
      WRITE(13,21)VS0,RS0,AAS,RC
*
      CALL NUCNAME(NINT(Z1),NINT(A1),ANAME1)
      CALL NUCNAME(NINT(Z2),NINT(A2),ANAME2)
      WRITE(13,*)'Nuclear System'
      WRITE(13,*) ANAME1,'  +  ', ANAME2	
      WRITE(13,*) ' '	
20    FORMAT(3(F8.3,4X),/)
21    FORMAT(4(F8.3,4X),/)
22    FORMAT(4(I2,5X),/)
*
      Z1Z2=Z1*Z2
      ARED=AVMASS*ACMASS/(AVMASS+ACMASS)/931.49386D0
*
      WRITE(13,*)'Radial integration parameters'
      WRITE(13,*) 'NPNTS      RMAX (fm)'
      WRITE(13,8) NPNTS,RMAX
8     FORMAT(I6,4X,F8.1,/)
*
* Setup nuclear + Coulomb + spin-orbit potential
*-----------------------------------------------
      CALL POTENT
*
* Integrate Schroedinger equation
*--------------------------------
      CALL BOUNDWAVE(N_0,L0,AJ0,ENER)
*
* calculate root-mean-square radius
*----------------------------------
      SUMR2=0.
      DO IR=1,NPNTS
         XR=IR*DR        
         SUMR2=SUMR2+XR**2*EIGWF(IR)**2
      ENDDO
      SUMR2=SQRT(SUMR2*DR)
*
* Output
*-------
      WRITE(6,*)' '
      WRITE(13,*)'# of nodes, orb ang mom, l+s ang mom, and energy'
      WRITE(6,*)'N_0     L0     AJ0     Energy  (MeV)'
      WRITE(13,*)'N_0     L0      AJ0     Energy  (MeV)'
      WRITE(6,104) N_0,L0,AJ0,ENER
      WRITE(13,104) N_0,L0,AJ0,ENER
      WRITE(13,*)'Wave function rms (fm) =',SUMR2
      WRITE(6,*)'Wave function rms (fm) =',SUMR2
 104  FORMAT(2(I2,6X),F4.1,6X,F10.4,/)	
      WRITE(13,*)'  '
      WRITE(13,*)'    R (fm)    x     Wavefunction u(r) = r R(r)'
      WRITE(13,*)' '
     
*
* output wavefunction
*--------------------
      DO IR=1,NPNTS/8.33,20
        XR=IR*DR        
        WRITE(13,*)XR,EIGWF(IR)
        WRITE(16,*)XR,EIGWF(IR)
      ENDDO
*
* just a reminder
*----------------
      PRINT*,' '
	PRINT*,'- Output of r x wavefunction in EIGEN.TXT and ',filena
	PRINT*,'- ',FILENA,'is ready for use by the other routines'
*
      RETURN
      END
*
*======================================================================
*
      SUBROUTINE BOUNDWAVE(NODE_0,L0,AJ,E)
*
* This subroutine solves the Schroedinger equation for bound states
*
      IMPLICIT REAL*8 (A-H,O-Z)
*
* NLG = maximum number of radial coordinate points
*------------------------------------------------- 
      PARAMETER(NLG =50000)
*
      COMMON /WAVE_FUNCTIONS/ EIGWF(NLG)
      COMMON /POTENTIAL/ U(NLG),WSO(NLG),UE(NLG)
     &                   ,VC(NLG)
      DIMENSION       WF0(NLG),
     &                VR(NLG),VE(NLG)
      COMMON /PAR/    DR,DR2,PI,FOUR_PI,NPNTS 
*
      NODE = NODE_0
      L    = L0
*
*  Find a good initial estimate of the eigenvalue
*------------------------------------------------
      IITER = 0
      NITER = 0
      INODE = NODE
      DEM    = 0.5D0
      COR0  = 0.5D0
      IF(L.EQ.0) COR0 =0.75D0
      COR   =   COR0+1.0D0/(DFLOAT(INODE)+1.0D0)/24.0D0
     *         -1.0D0/(DFLOAT(INODE)+1.0D0)**2/48.0D0
      A12   =    1.0D0/12.0D0
      ACCUR = 0.1D0
      ACCUR0= 0.1D-8
*
* Find bottom of the well
*------------------------
      AL  = DFLOAT( L*(L+1) )
      ALS = AJ*(AJ+1.0D0)-AL-0.75D0
      EMIN = 0.0D0
      DO IR=1,NPNTS
        XR = IR
        VE(IR) = UE(IR)
        VR(IR) = U(IR) + AL/IR**2 + ALS*WSO(IR)
        IF( EMIN .GT. ( VR(IR)/VE(IR) ) ) EMIN = VR(IR)/VE(IR)
      ENDDO
      IF( NITER. GT. 6 ) GO TO 100
      IF( E .LE. EMIN  .OR.  E .GE. 0.0D0 ) E = EMIN*0.5D0
*
* Find turning point and WKB eigen-energy
*----------------------------------------
9     NITER=NITER+1
      S=0.0D0
      S1=0.0D0
      DO 1 IR=1,NPNTS
        P2 = E*VE(IR)-VR(IR)
        IF ( P2.LE. 0.0D0) GO TO 1
        P  = DSQRT( P2 )
        S  = S+P
        S1 = S1+VE(IR)/P
        ITURN = IR
1     CONTINUE
      S = S-PI*( DFLOAT(INODE)+COR )
      DE=-S/S1*2.0D0
108   IF( (E+DE) .GT. EMIN .AND. (E+DE) .LT. 0.0D0 )  GO TO 109
      DE=DE*0.5D0
      GO TO 108
109   E=E+DE
      IF( NITER .GT. 50 ) GO TO 1000
      IF( DABS(DE) .GT. ACCUR ) GO TO 9
      XR=ITURN*DR
100   IITER=IITER+1
110   CONTINUE
*
*  Integrate outward from origin to the turning point
*----------------------------------------------------
      A     = 0.0D0
      B     = 1.0D-10
      WF0(1) = B
      D     = VR(1)-VE(1)*E
      B     = B*(1.0D0-A12*D)
      SEI   = WF0(1)**2*VE(1)
      ICHECK=0
      DO 101 IR=2,ITURN
        C   = 2.0D0*B-A+WF0(IR-1)*D
        A   = B
        B   = C
        D   = VR(IR)-VE(IR)*E
        C   = C/(1.0D0-A12*D)
        SEI = SEI+C**2*VE(IR)
        WF0(IR) = C
        IF(IR.LT.5)  GO TO 101
        IF( WF0(IR)*WF0(IR-1) .LT. 0.0D0 .OR. WF0(IR) .EQ. 0.0D0)
     &   ICHECK=ICHECK+1
101   CONTINUE
      IF(ICHECK.GT.INODE) THEN
        E=E*1.2D0
        DO 1200 IR=1,NPNTS
          IF( VE(IR)*E .GT. VR(IR) ) ITURN=IR
1200    CONTINUE
        GO TO 110
      ENDIF
      IF(ICHECK.LT.INODE) THEN
        E=E*0.9D0
        DO 1201 IR=1,NPNTS
          IF( VE(IR)*E .GT. VR(IR)) ITURN=IR
1201    CONTINUE
        GO TO 110
      ENDIF
      SEI = SEI/C**2-0.5D0*VE(ITURN)
      WFI = C
      WF1I=(274.0D0*WF0(ITURN)-600.0D0*WF0(ITURN-1)
     &  +600.0D0*WF0(ITURN-2)
     &    -400.0D0*WF0(ITURN-3)+150.0D0*WF0(ITURN-4)
     &      -24.0D0*WF0(ITURN-5))/(120.0D0*DR)
      DWFI = WF1I/WFI
      SEI = SEI/DR+A12*DWFI*DR**2
*
* Integrate inward from infinity to the turning point
*----------------------------------------------------
      D = VR(NPNTS)-VE(NPNTS)*E
      A = DEXP(-DFLOAT(NPNTS)*DSQRT(DABS(D)))
      IF ( A .LT. 1.0D-12) A =0.0D0
      WF0(NPNTS) = A
      A = A*(1.0D0-A12*D)
      D = VR(NPNTS-1)-VE(NPNTS-1)*E
      B = DEXP(-DFLOAT(NPNTS-1)*DSQRT(DABS(D)))
      IF ( B .LE. 1.0D-12 ) B = 1.0D-12
      WF0(NPNTS-1) = B
      SEO = B**2*VE(NPNTS-1)
      B = B*(1.0D0-A12*D)
      DO 201 IR=NPNTS-2,ITURN,-1
        C   = 2.0D0*B-A+WF0(IR+1)*D
        A   = B
        B   = C
        D   = VR(IR)-VE(IR)*E
        C   = C/(1.0D0-A12*D)
        SEO = SEO+C**2*VE(IR)
        WF0(IR)=C
201   CONTINUE
      WFO = C
      WF1O=(-274.0D0*WF0(ITURN)+600.0D0*WF0(ITURN+1)
     *  -600.0D0*WF0(ITURN+2)
     *   +400.0D0*WF0(ITURN+3)-150.0D0*WF0(ITURN+4)
     *     +24.0D0*WF0(ITURN+5))/(120.0D0*DR)
      DWFO = WF1O/WFO
      SEO = SEO/WFO**2-0.50D0*VE(ITURN)
      SEO = SEO/DR-A12*DWFO*DR**2
      SE  = DWFI-DWFO
*
* Compute the correction to the eigenvalue
*-----------------------------------------
      DE = SE/(SEI+SEO)
221   IF( DABS(DE) .LT. DEM ) GO TO 222
      DE = DE*0.8D0
      GO TO 221
222   E1 = E+DE
      IF( E1 .GE. 0.0D0 .OR. E1 .LE. EMIN ) THEN
        DE = DE*0.9D0
        GO TO 222
      ENDIF
      E = E1
      IF( IITER .GT. 200 ) GO TO 1000
      XTEST = DABS(DE)/(DABS(E)+1.D-25)
      IF( XTEST .GT. ACCUR0 ) GO TO 100
*
* Match the wavefunction
*-----------------------
      SE = WFI/WFO
      DO 303 IR=ITURN,NPNTS
303   WF0(IR)=WF0(IR)*SE
*
* Normalize the wavefunction
*---------------------------
      SN = 0.0D0
      DX = DR/15.0D0
      AX = DX*7.0D0
      NUP= NPNTS
      IF( MOD(NPNTS,2) .NE.0 ) NUP=NPNTS-1
      DO 301 IR=NUP,1,-1
        SN = SN+WF0(IR)**2*VE(IR)*AX
        AX = DR+DX
        DX =-DX
301   CONTINUE
      SN=1.0D0/DSQRT(SN)
      DO 302 IR=1,NPNTS
        AX = WF0(IR)*SN*DSQRT(VE(IR))
        XR = IR*DR
      WF0(IR) = AX/XR
302   EIGWF(IR) = AX
      RETURN
*
* OOPS! Bad luck happens.
*------------------------
1000  WRITE(6,*) 'The iteration procces failed  NITER=',NITER,'
     &    IITER='
     *     ,IITER,'  NODE=',NODE,'  E=',E
*
      STOP
      END
*
*======================================================================
*
      SUBROUTINE POTENT
*
*  Construct potentials for bound state wavefunctions
*
      IMPLICIT REAL*8 (A-H,O-Z)
*
* NLG = maximum number of radial coordinate points
*------------------------------------------------- 
      PARAMETER(NLG = 50000)
*
      COMMON /POTENTIAL/ U(NLG),WSO(NLG),UE(NLG)
     &                   ,VC(NLG)
      COMMON /PAR/     DR,DR2,PI,FOUR_PI,NPNTS 
      COMMON /POT/ V0,VS0,R0,RS0,RC,AA,AAS,Z1Z2,ARED
*
      NCOUNT=0
*
* Construct the Coulomb potential
*--------------------------------
      DR2 = DR**2
      B_N = 931.49386D0*ARED	
      B_MASS= (197.327D0)**2/(2.0D0*B_N)	
      CHARGE = 1.43996D0*Z1Z2				
      DO IR=1,NPNTS
        XR = IR*DR
        XRI=(XR-R0)/AA
        IF(XRI.GT.30.0D0)XRI=30.0D0
        EXT=DEXP( XRI )
        XRIS=(XR-RS0)/AAS
        IF(XRIS.GT.30.0D0)XRIS=30.0D0
        EXTS=DEXP( XRIS )
*
* Central
*--------
        U(IR) = V0/( 1.0D0+EXT )/B_MASS*DR2
*
* Spin-orbit
*----------- 
        RCOMP=1.41382     ! pion Compton wavelength
        WSO(IR) = VS0*RCOMP**2*EXTS/((1.0D0+EXTS)**2*XR*AAS)
     &                  /B_MASS*DR2/2.
*
* Scaling parameter
*------------------
        UE(IR) = DR2/B_MASS
*
* Coulomb potential of a uniformly distributed charge
*----------------------------------------------------
        VC(IR)=CHARGE/XR
        IF(XR.LT.RC) VC(IR)=CHARGE/RC*
     &      (1.5D0-0.5D0*(XR/RC)**2)
        VC(IR)=VC(IR)*DR2/B_MASS
*
* add to WS potential
*--------------------
        U(IR)=U(IR)+VC(IR)
      ENDDO
*
      RETURN
      END
*
*====================================================================
*
      SUBROUTINE  SMAT(AP,ZP,AT,ZT,ECA)
*
* Computes the eikonal S-matrices for heavy ion collisions.
* The S-matrices are functions of the impact parameter B. 
* The program is appropriate for heavy ion bombarding energies 
* E_lab > 30 MeV/nucleon.
* 
     	IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NBIG=10000,NRMAX=300,NB=300,NQ=9999)
*
     	COMPLEX*16 CI,WN(NBIG),
     &           SMATRX(NBIG),PHASEN(NBIG),PSHIFT
      REAL*8 RR(NBIG),DENSP(NBIG),DENST(NBIG),REALWN(NBIG),DIMAWN(NBIG),
     &       Q(NBIG),B(NBIG),QDENSP(NBIG),QDENST(NBIG),
     &       FELAST(NBIG)
*
      CHARACTER*12 FILENA,FILENA2
      CHARACTER NAMEP*5,NAMET*5
*
      COMMON/A1/ BETA,GAMMA
      COMMON /CONST1/ PI,FOUR_PI,HBARC,AMU,FINE,E2
      COMMON /CONST2/ CI
      COMMON /CUT/ BMAX,RMAX
*
      OPEN(UNIT=17,FILE='SMAT.TXT', STATUS='UNKNOWN')
*
      PRINT*,' Enter name of OUTPUT file'
      READ(*,*)FILENA
*
      OPEN(UNIT=14,FILE=FILENA, STATUS='UNKNOWN')
*
* read option: IOPT=0 for t-rho-rho, 1 for input opt. potential
*--------------------------------------------------------------
      CALL SKIPCOM(12)
      READ(12,*)IPOT
      NT=NINT(AT-ZT)
      NP=NINT(AP-ZP)
*
* mesh in momentum, impact parameter, and angle
*----------------------------------------------
      QMAX  =  6.   ! maximum value of momentum transfer
      DQ=QMAX/NQ
     	DO I=1,NQ
        Q(I)=I*DQ
      ENDDO
*
* resize R and B mesh
*--------------------		      
      DR=RMAX/NRMAX
      DO I=1,NRMAX
        RR(I)=I*DR
      ENDDO
*
     	DB=BMAX/NB
     	DO I=1,NB
        B(I)=I*DB
      ENDDO
*
* energies and momenta of the nuclei
*-----------------------------------
      ELAB=ECA*AP
      AMASS1=EXCESS(ZP,AP,1.D0) + AP*AMU - ZP*0.510999
      AMASS2=EXCESS(ZT,AT,1.D0) + AT*AMU - ZT*0.510999
      REDM=AMASS1*AMASS2/(AMASS1+AMASS2)/AMU
     	ECM=ELAB*REDM/AP
     	ETA=ZT*ZP*1.44/(HBARC*BETA)
*
* output of reaction data
*------------------------
     	WRITE(17,99)
      CALL NUCNAME(NINT(ZP),NINT(AP),NAMEP)
      CALL NUCNAME(NINT(ZT),NINT(AT),NAMET)
      WRITE(17,113)NAMEP,'+',NAMET
      WRITE(17,114)ECA	
*
* I/O format
*-----------
  100	FORMAT(30A1)
  99	FORMAT(2X,'**** Output of routine SMAT *****',/,
     &       2x,'*     (units in fm and MeV)     *',/,
     &       2x,'*********************************',/) 
  113	FORMAT(2X,5A,5A) 
  114	FORMAT(2X,'Elab (MeV/nuc)=',f9.4,/) 
  104	FORMAT(2X,'Nucleon-nucleon scattering data:')
  105	FORMAT(2X,'isospin averaged NN-cross section (mb) =',f8.4,/,
     &      2x,'imaginary coefficient of the f_NN amplitude =',f8.4,/,
     &      2x,'exponential coefficient of the f_NN amplitude =',f8.4,/)
*
* read parameters and build up Woods-Saxon optical potential
*-----------------------------------------------------------
      IF(IPOT.EQ.0) THEN
        CALL SKIPCOM(12)
        READ(12,*)V0_WS,R_WS,D_WS,VI_WS,RI_WS,DI_WS
        CALL OMP_WS(RR,NRMAX,V0_WS,VI_WS,R_WS,RI_WS,D_WS,DI_WS,WN)
      ENDIF
* or read optical potential from files and interpolate to NRMAX points
*---------------------------------------------------------------------
      IF(IPOT.EQ.1)THEN
        CALL SKIPCOM(12)
        READ(12,*)FILENA,FILENA2
        CALL FINFOUT(FILENA,NRMAX,DR,REALWN)
        CALL FINFOUT(FILENA2,NRMAX,DR,DIMAWN)
        DO I=1,NRMAX
          WN(I)=REALWN(I)+CI*DIMAWN(I)
        ENDDO
      ENDIF
*
* phase-shifts from input optical potentials 
*-------------------------------------------
      IF(IPOT.LE.1) THEN
        CALL PHNUC(RR,NRMAX,WN,B,NB,PHASEN)
        DO I=1,NRMAX
          SMATRX(I)=EXP(CI*PHASEN(I))
        ENDDO
      ENDIF
*    
*  nucleon-nucleon cross section (in fm^2)
*  with Pauli blocking correction (IPAULI=1)
*
*  For no correction, set IPAULI=0
*---------------------------------
      CALL SKIPCOM(12)
      READ(12,*)IPAULI
      CALL SIGNNE(AP,ZP,AT,ZT,ECA,IPAULI,SIGNN)
      SIGNN=SIGNN/10.
*
*  real part of the nucleon-nucleon scattering 
*  amplitude and momentum dependence 
*-----------------------------------
      CALL PHNNE(AP,ZP,AT,ZT,ECA,BNN,ALFANN)
*
      WRITE(17,104)
      WRITE(17,105)SIGNN*10.,BNN,ALFANN
*
* read density parameters (ALPHA is a correction for 
* the proton size)
*
* For input, see H. de Vries et al., Atom. Dat. Nucl. Dat. Tables,
* Vol 36, 495 (1987)
*-------------------
      CALL SKIPCOM(12)
      READ(12,*)IDPROJ,IDTARG
* projectile density
*-------------------
      IF(IDPROJ.EQ.1)  CALL GAUSS(NRMAX,DR,DENSP)  
      IF(IDPROJ.EQ.2)  CALL YUKAWA(NRMAX,DR,DENSP)
      IF(IDPROJ.EQ.3)  CALL FERMI(NRMAX,DR,DENSP)
      IF(IDPROJ.EQ.4)  CALL SAXON(NRMAX,DR,DENSP)
      IF(IDPROJ.EQ.5)  CALL DROP(AP,DFLOAT(NP),NRMAX,DR,DENSP)
      IF(IDPROJ.EQ.10) THEN
        CALL SKIPCOM(12)
        READ(12,*)FILENA,ALPHAP
        CALL FINFOUT(FILENA,NRMAX,DR,DENSP)
      ENDIF
* target density
*---------------
      IF(IDTARG.EQ.1)  CALL GAUSS(NRMAX,DR,DENST)  
      IF(IDTARG.EQ.2)  CALL YUKAWA(NRMAX,DR,DENST)
      IF(IDTARG.EQ.3)  CALL FERMI(NRMAX,DR,DENST)
      IF(IDTARG.EQ.4)  CALL SAXON(NRMAX,DR,DENST)
      IF(IDTARG.EQ.5)  CALL DROP(AT,DFLOAT(NT),NRMAX,DR,DENST)
      IF(IDTARG.EQ.10) THEN
        CALL SKIPCOM(12)
        READ(12,*)FILENA,ALPHAT
        CALL FINFOUT(FILENA,NRMAX,DR,DENST)
      ENDIF
*
* Fourier transform of the densities
*-----------------------------------
      CALL FOURIER0(DENSP,RR,NRMAX,QDENSP,Q,NQ)     
      CALL FOURIER0(DENST,RR,NRMAX,QDENST,Q,NQ)
* fold with proton densities
*---------------------------
      IF(IDPROJ.EQ.10) THEN
        DO IQ=1,NQ
          QDENSP(IQ)=QDENSP(IQ)*EXP(-Q(IQ)**2*ALPHAP**2/4.)
        ENDDO
      ENDIF
      IF(IDTARG.EQ.10) THEN
        DO IQ=1,NQ
          QDENST(IQ)=QDENST(IQ)*EXP(-Q(IQ)**2*ALPHAT**2/4.)
        ENDDO
      ENDIF
* normalize distributions to nucleon number
*------------------------------------------
      APDUM=QDENSP(1)
      ATDUM=QDENST(1)
      DO IQ=1,NQ
        QDENSP(IQ)=AP/APDUM*QDENSP(IQ)
        QDENST(IQ)=AT/ATDUM*QDENST(IQ)
      ENDDO       
*
*  Eikonal phase-shifts and S-matrix
*-----------------------------------
     	FNR=ALFANN*SIGNN/(4.*PI)
     	FNI=SIGNN/(4.*PI)
*
* this part is for phase-shifts from nucleon-nucleon collisions
* integrate over momentum Q for a set of impact parameters B 
*-----------------------------------------------------------
     	DO IB=1,NB
* 
* integration over 'q' - Simpson's rule
        DO IQ=1,NQ
          QB=Q(IQ)*B(IB)
          FJ0=BESSJ0(QB)
          EF=EXP(-BNN*Q(IQ)**2)
          FELAST(IQ)=Q(IQ)*QDENSP(IQ)*QDENST(IQ)*FJ0*EF
        ENDDO
        CALL RSIMP(FELAST,DUM,DQ,NQ)
* end of integral over 'q'
        PSHIFT=DUM*(FNR+CI*FNI)
        SMATRX(IB)=EXP(CI*PSHIFT)
      ENDDO
*
* print out S-matrices
*---------------------
  11 	DO JB=1,NB
        WRITE(17,3000)B(JB),DREAL(SMATRX(JB)),DIMAG(SMATRX(JB))
        WRITE(14,3000)B(JB),DREAL(SMATRX(JB)),DIMAG(SMATRX(JB))
      ENDDO        
*
 2000 FORMAT(2X,2(G10.4,14X))   
 3000 FORMAT(2X,3(G10.4,2X))   
*
* just a reminder
*----------------
      PRINT*,' '
	PRINT*,'- Output of S-Matrix in SMAT.TXT and ',FILENA
	PRINT*,'- ',FILENA,' is ready for use by the other routines.'
*
 1000 STOP
      END
*
*====================================================================
*
      SUBROUTINE OMP_WS(R,NR,V0_WS,VI_WS,R_WS,RI_WS,D_WS,DI_WS,WN)
*
* Woods-Saxon optical model potential. 
*
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 R(NR)
      COMPLEX*16 WN(NR),CI
      COMMON/CONST2/CI
*
      DO I=1,NR
        VDUM1 = -V0_WS/(1.+DEXP((R(I)-R_WS)/D_WS))
        VDUM2 = -VI_WS/(1.+DEXP((R(I)-RI_WS)/DI_WS))
        WN(I) = VDUM1+CI*VDUM2
      END DO     
* print potential, if requested
*        WRITE(7,'(A12,2X,(A12,1X,A12))')
*     &    " r[fm]"," Re{U} [MeV]"," Im{U} [MeV]"
*        DO I=0,NR
*          WRITE(7,'(F12.4,2X,(E12.4,1X,E12.5))')R(I),WN(I)
*        END DO
*
      RETURN 
      END
*
*====================================================================
*
      SUBROUTINE CROSS(L,M,IOPT,NR,NB,NK)
*
* If IOPT = 1
* Output is d sigma / d p_z 
*
* If IOPT = 2
* Output is d sigma / d p_t^2 
*
* If IOPT = 3
* Output is d sigma / d p_y 
*
* If IOPT = 4
* Output is d sigma / d p_z d p_t, where p_z = IKZ * DK and 
* p_t = IKT * DK are the longitudinal and transverse momentum 
*
* If IOPT = 5
* Output is core+target elastic scattering differential cross section 
*
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NLG=310,NWF=300,NGAUSS=40,NDAT=300)
      REAL*8 FAK(0:40),FINT(0:NR),SIG(0:NK),AUX1(NWF),RW(NWF),
     &   R(0:NR),KZ(0:NK),KT(0:NK),SIMPB(0:NB),IV(0:40),BW(NWF),
     &   BN(0:NB),WF(0:NR),SN(0:NB),SIMPR(0:NR),AUX2(NWF),
     &   DUM1(0:NK,0:NK),BLMSQ(0:NK,0:NB),SIMPK(0:NK),
     &   BETAN(0:NGAUSS),SINT(0:NB,0:NR),SC(0:NB),
     &   SCR(0:NB),SCI(0:NB),SNR(0:NB),SNI(0:NB),
     &   SIGEL(0:NK),AUX3(0:NB)
      COMPLEX*16	CI,ZINT(0:NK,0:NR),SUMEL,ALPHAC(NGAUSS),FC,TEMP
      COMMON /GF1 / IV
      COMMON /GF2 / FAK
      COMMON /BT/ BETAN
      COMMON /GS1/ NGS
      COMMON /GS2/ ALPHAC
      COMMON /COUL/ PCM,PCM2,ETA,ETA2
      COMMON /PMOM/ PMAX
      COMMON /CONST1/ PI,FOUR_PI,HBARC,AMU,FINE,E2
      COMMON /CONST2/ CI
      COMMON /SYST/ APROJ,ZPROJ,ATARG,ZTARG,ACORE,ZCORE,ECA,S2C
      COMMON /ML/ MLML
      COMMON /MLC/ CLM
      COMMON /CUT/ BMAX,RMAX
*
      CHARACTER*12 FILENA,FILEWS,FILESC,FILESN
*
* input files
*------------
      PRINT*,' Enter name of wavefunction input file'
      READ(*,*)FILEWS
      PRINT*,' Enter name of S-matrix (valence nucleon+target) file'
      READ(*,*)FILESN
      PRINT*,' Enter name of S-matrix (core+target) file'
      READ(*,*)FILESC
      PRINT*,' Enter name of output file'
      READ(*,*)FILENA
*
      OPEN(UNIT=20,FILE=FILENA,STATUS='UNKNOWN')
      OPEN(UNIT=21,FILE=FILEWS,STATUS='UNKNOWN')
      OPEN(UNIT=22,FILE=FILESC,STATUS='UNKNOWN')
      OPEN(UNIT=23,FILE=FILESN,STATUS='UNKNOWN')
*
      TOL=1.D-20
      M=IABS(M)
*
* Parameters for the expansion in Gaussians (in fm)
*--------------------------------------------------
      BETA0=50.
      NGS=NGAUSS
      IF(ETA.EQ.0.) THEN
     	  BETA0=20.
        NGS=NGAUSS/2
      ENDIF				 
      BETAN(0)=1.D10
      DO  I=1,NGS
        BETAN(I)=BETA0/I
      ENDDO
*
* read Gaussian expansion coefficients for S-matrix (core+target)
*----------------------------------------------------------------
      CALL SMAT_FIT(IOPT,NGS,BETA0,NDAT,ALPHAC)
*
* read radial wavefunction
*-------------------------
      SUM=0.
      DO I=1,NWF
        READ(21,*)RW(I),AUX1(I)
        AUX1(I)=AUX1(I)/RW(I)   !change u(r) ---> R(r)
        SUM=SUM+RW(I)**2*AUX1(I)**2 
      ENDDO
      DROLD=RW(4)-RW(3)
      SUM=SUM*DROLD
* check norm of the wavefunction
*-------------------------------
      IF(ABS(SUM-1.).GT.0.05)THEN
      	PRINT*,'STOP: norm of the wavefunction is not correct'
          PRINT*,'NORM = ', SUM 
          STOP
      ENDIF	
      CLOSE(UNIT=52)
* rescale mesh in radial coordinate
*----------------------------------
      DR=RMAX/NR
      R(0)=1.D-10
      DO I=1,NR
        R(I)=I*DR
      ENDDO
*
* interpolate wavefunction in the new radial mesh
*------------------------------------------------
      DO I=1,NR
        WF(I)=BETWEEN(R(I),AUX1,RW,NWF,NWF)
      ENDDO
      WF(0)=0.
*
* read S-matrices for nucleon+target
*-----------------------------------
      DO I=1,NWF
       READ(23,*)BW(I),AUX1(I),AUX2(I)
      ENDDO
*
* rescale mesh in impact parameter
*---------------------------------
      DB=BMAX/NB
      DO I=0,NB
        BN(I)=I*DB
      ENDDO
*
* interpolate S_n in the new mesh
*--------------------------------
      DO IB=1,NB
        SNR(IB)=BETWEEN(BN(IB),AUX1,BW,NWF,NWF)
        SNI(IB)=BETWEEN(BN(IB),AUX2,BW,NWF,NWF)
      ENDDO
      SNR(0)=0.
      SNI(0)=0.
*
* add Coulomb phase
*------------------
      DO IB=1,NB
        SN(IB)=SQRT(SNR(IB)**2+SNI(IB)**2)
        SCOULR=COS(2.*ETA2*DLOG(PCM2*BN(IB)))
        SCOULI=SIN(2.*ETA2*DLOG(PCM2*BN(IB)))
        SNRR=SNR(IB) 
        SNII=SNI(IB) 
        SNR(IB)=SNRR*SCOULR-SNII*SCOULI        
        SNI(IB)=SNII*SCOULR+SNRR*SCOULI        
      ENDDO
*
* read (interpolate) S-matrix for core+target
*--------------------------------------------
      OPEN(UNIT=22,FILE=FILESC,STATUS='UNKNOWN')
      DO I=1,NWF
       READ(22,*)BW(I),AUX1(I),AUX2(I)
      ENDDO
* interpolate S_c in the new B-mesh
*----------------------------------
      DO IB=1,NB
        SCR(IB)=BETWEEN(BN(IB),AUX1,BW,NWF,NWF)
        SCI(IB)=BETWEEN(BN(IB),AUX2,BW,NWF,NWF)
      ENDDO
      SCR(0)=0.
      SCI(0)=0.
* add Coulomb phase
*------------------
      DO IB=1,NB
        SC(IB)=SQRT(SCR(IB)**2+SCI(IB)**2)
        SCOULR=COS(2.*ETA*DLOG(PCM*BN(IB)))
        SCOULI=SIN(2.*ETA*DLOG(PCM*BN(IB)))
        SCRR=SCR(IB)
        SCII=SCI(IB)
        SCR(IB)=SCRR*SCOULR-SCII*SCOULI        
        SCI(IB)=SCII*SCOULR+SCRR*SCOULI        
      ENDDO
*
* mesh in momentum
*-----------------
      PMAX=PMAX/HBARC
      DK=PMAX/NK
      DO I=0,NK
        KZ(I)=I*DK
        KT(I)=I*DK
      ENDDO
*
* Simpson coefficients for integral over BN
*------------------------------------------
      SIMPB(0)=1./3.*DB
      SIMPB(NB)=1./3.*DB
      IS=4
      DO K=1,NB-1
        SIMPB(K)=IS/3.*DB
        IS=6-IS
      ENDDO
*
* Simpson coefficients for integrals over momentum
*-------------------------------------------------
      SIMPK(0)=1./3.*DK
      SIMPK(NK)=1./3.*DK
      IS=4
      DO K=1,NK-1
        SIMPK(K)=IS/3.*DK
        IS=6-IS
      ENDDO
*
* Simpson coefficients for z- and rho-integration
*------------------------------------------------
      SIMPR(0)=1./3.*DR
      SIMPR(NR)=1./3.*DR
      IS=4
      DO K=1,NR-1
        SIMPR(K)=IS/3.*DR
        IS=6-IS
      ENDDO
*
* do integral along the z-direction and store in FINT
*----------------------------------------------------
      CALL FINTEG(L,M,R,NR,WF,SIMPR,FINT)
*
* do integral of |S_c|^2 over PHI and store it in SINT
*-----------------------------------------------------
      CALL PINTEG(R,NR,BN,NB,SC,SINT)
*
* integrate over BN and RHO
*--------------------------
      STRIP=0. 
      DO IB=1,NB
        SUMRHO=0.
        DO IR=0,NR
          RHO=R(IR)
          SUMRHO=SUMRHO+SIMPR(IR)*RHO*FINT(IR)*SINT(IB,IR)
        ENDDO
        STRIP=STRIP+SIMPB(IB)*BN(IB)*(1.-SN(IB)**2)*CLM**2*
     &     	  2.*PI*SUMRHO/(2.*L+1)      
      ENDDO
*
* add diffraction dissociation X-section
*---------------------------------------
      CALL COMPDIF(L,M,R,NR,BN,NB,SNR,SNI,SCR,SCI,WF,SIMPR,AUX3)
      SUMDIF=0.
      DO IB=0,NB
        SUMDIF=SUMDIF+SIMPB(IB)*BN(IB)*AUX3(IB)
      ENDDO
      SDIF=2.*PI*SUMDIF/(2.*L+1.)
      STOT=S2C*(STRIP+SDIF)   
      STRIP=MLML*STRIP 	  
      SDIF=MLML*SDIF 	  
*
* Calculate longitudinal momentum distribution
*=============================================
      IF(IOPT.EQ.1) THEN
* do integral along the z-direction and store in ZINT
*----------------------------------------------------
        CALL ZINTEG(L,M,R,NR,WF,KZ,NK,SIMPR,ZINT)
*
* momentum distribution matrix
*-----------------------------
        CALL COMPL(R,NR,NK,NB,SIMPR,ZINT,SINT,BLMSQ)
*
* integral over KZ (from 0 to infinity and multiply by 2)
*--------------------------------------------------------
        SUMTOT=0. 
        DO IKZ=0,NK
* integral over BN
*-----------------
          SUM=0.
          DO IB=1,NB
            SUM=SUM+SIMPB(IB)*BN(IB)*(1.-SN(IB)**2)*CLM**2*
     &     		  BLMSQ(IKZ,IB)      
          ENDDO
          SIG(IKZ)=2.*PI/(2.*L+1.)*SUM          
          SUMTOT=SUMTOT+2.*SIMPK(IKZ)*SIG(IKZ)			
        ENDDO
* loop over KZ
*-------------
        DO IKZ=NK,0,-1
          DUM=SIG(IKZ)*STOT/SUMTOT
          WRITE(20,100)-KZ(IKZ)*HBARC,10*MLML*DUM/HBARC
        ENDDO
* use reflection symmetry around KZ=0
*------------------------------------
        DO IKZ=1,NK
          DUM=SIG(IKZ)*STOT/SUMTOT
          WRITE(20,100)KZ(IKZ)*HBARC,10*MLML*DUM/HBARC
        ENDDO
* output of total cross section
*------------------------------
        STOT=MLML*STOT 
        PRINT*,' '
        WRITE(*,111)10.*STRIP,10.*SDIF,10.*STOT	          
        PRINT*,' '
        PRINT*,' Longitud. momentum distrib. ds/dp_z (mb/MeV/c) in ',
     &	  FILENA
      ENDIF
111   FORMAT(' Stripping X-section (M=-L + M=L) [mb]  =  ',G10.4,/,
     & ' Diffraction X-section (M=-L + M=L) [mb]  =  ',G10.4,/,
     & ' Total knockout X-section (M=-L + M=L) [mb]  =  ',G10.4)
*
* Calculate transverse momentum distribution
*===========================================
*
      IF(IOPT.EQ.2.OR.IOPT.EQ.3) THEN
* sum over Gaussians
*-------------------
        CALL GAUST(L,M,R,NR,KT,NK,BN,NB,NGS,SIMPR,WF,BLMSQ)
*
* Calculate d sigma / d k_t^2
*---------------------------- 
        IF(IOPT.EQ.2) THEN
* loop and integral over KT
*--------------------------
          SUMTOT=0. 
          DO IKT=0,NK
* integral over BN
*-----------------
            SUM=0.
            DO IB=0,NB
              SUM=SUM+SIMPB(IB)*BN(IB)*(1.-SN(IB)**2)*CLM**2*
     &     		  BLMSQ(IKT,IB)    
            ENDDO
            SIG(IKT)=2.*PI/(2.*L+1.)*SUM
            SUMTOT=SUMTOT+SIMPK(IKT)*KT(IKT)*SIG(IKT)			
          ENDDO
          SUMTOT=2.*PI*SUMTOT
* loop over KT
*-------------
          DO IKT=0,NK
            DUM=SIG(IKT)*STOT/SUMTOT
            IF(DUM.EQ.0.)DUM=1.D-30
            WRITE(20,100)KT(IKT)*HBARC,10*MLML*DUM/HBARC/HBARC		       		 
          ENDDO
* output of total cross section
*------------------------------
          STOT=MLML*STOT 
          PRINT*,' '
          WRITE(*,111)10.*STRIP,10.*SDIF,10.*STOT	          
          PRINT*,' '
          PRINT*,' Transv. mom. dist. ds/d p_t^2 (mb/(MeV/c)^2) in '
     & ,FILENA
        ENDIF
*
* Calculate d sigma / d k_y
*-------------------------- 
        IF(IOPT.EQ.3) THEN
*
* loop and integral over KY
*--------------------------
          SUMY=0.
          DO IKY=0,NK
            PKY=KT(IKY)
* loop and integral over KX
*--------------------------
            SUMX=0.
            DO IKX=0,NK
              PKX=KT(IKX)
              PKT=SQRT(PKX*PKX+PKY*PKY)
              KMID=INT(PKT/DK)	
              AMID=0.
              AUX=0.										              		          
              IF(KMID.LT.NK)AMID=1.-(PKT-KT(KMID))/DK
* integral over BN
*-----------------
              SUMB=0.
              DO  IB=0,NB
                IF(KMID+1.LE.NK.AND.KMID.GE.0) THEN
                  AUX=AMID*BLMSQ(KMID,IB)+(1.-AMID)*BLMSQ(KMID+1,IB)
                  SUMB=SUMB+SIMPB(IB)*BN(IB)*(1.-SN(IB)**2)*CLM**2*
     &               AUX
                ENDIF 
              ENDDO     
              SEC=2.*PI/(2.*L+1.)*SUMB
              SUMX=SUMX+SIMPK(IKX)*SEC
            ENDDO	
            SIG(IKY)=SUMX		  	            		  							 
            SUMY=SUMY+SIMPK(IKY)*SIG(IKY)
          ENDDO
          SUMY=4.*SUMY
* loop over KY
*-------------
          DO IKY=NK,0,-1
            SIG(IKY)=2.*SIG(IKY)*STOT/SUMY
            IF(SIG(IKY).EQ.0.)SIG(IKY)=1.D-30
            WRITE(20,100)-KT(IKY)*HBARC,10.*MLML*SIG(IKY)/HBARC
          ENDDO
* use reflection symmetry around KY=0
*------------------------------------
          DO IKY=1,NK
            WRITE(20,100)KT(IKY)*HBARC,10*MLML*SIG(IKY)/HBARC
          ENDDO
          STOT=MLML*STOT 
          PRINT*,' '
          WRITE(*,111)10.*STRIP,10.*SDIF,10.*STOT	          
          PRINT*,' '
          PRINT*,' Transv. mom. dist. ds/d p_y (mb/MeV/c) in ',
     &     		FILENA
        ENDIF				  		  		       		 
      ENDIF
*
* Calculate double-differential momentum distribution
*====================================================
*
      IF(IOPT.EQ.4) THEN         
*
* do integral along the z-direction and store in ZINT
*----------------------------------------------------
        CALL ZINTEG(L,M,R,NR,WF,KZ,NK,SIMPR,ZINT)
* sum over Gaussians
*-------------------
        CALL 
     &  GAUSLT(M,R,NR,KT,NK,BN,NB,SN,NGS,SIMPR,SIMPB,ZINT,DUM1)
* integral over KZ
*-----------------
        SUMZ=0.
        DO IKZ=0,NK
          SUMT=0.
* integral over KT
*-----------------
          DO IKT=0,NK
            DUM1(IKZ,IKT)=CLM**2*2*PI*KT(IKT)/(2*L+1.)*DUM1(IKZ,IKT)
            SUMT=SUMT+SIMPK(IKT)*DUM1(IKZ,IKT)
          ENDDO
          SUMZ=SUMZ+2.*SIMPK(IKZ)*SUMT
        ENDDO
        DO IKZ=0,NK
          DO IKT=0,NK
            DUM1(IKZ,IKT)=DUM1(IKZ,IKT)*STOT/SUMZ
            IF(DUM1(IKZ,IKT).EQ.0.)DUM1(IKZ,IKT)=1.D-30
            WRITE(20,90)KZ(IKZ)*HBARC,KT(IKT)*HBARC,
     &	    		 10*MLML*DUM1(IKZ,IKT)/(HBARC)**2
          ENDDO
        ENDDO
        STOT=MLML*STOT 
        PRINT*,' '
        WRITE(*,111)10.*STRIP,10.*SDIF,10.*STOT	          
        PRINT*,' '
        PRINT*,' Double dif. X-section ds/dp_z dp_t (mb/(MeV/c)^2) in'
     & , FILENA
      ENDIF
*
* Calculate elastic scattering cross section of core+target d sigma/ d p_y
*=========================================================================
      IF(IOPT.EQ.5) THEN   
*
* zeroth Coulomb phase
*---------------------
        CALL CP0(ETA,ARGAM)
* 
        NBIG=1000  
        DBSM=BMAX/NBIG
        SUMTOT=0.
        IG=4
        DO IKT=0,NK
          SUMEL=0.
* impact parameter integral with a fine mesh
*-------------------------------------------
          DO IB=0,NBIG
            BNEW=IB*DBSM		 
            IOLD=BNEW/DB
            TEMP=SCR(IOLD)-CI*SCI(IOLD)            
            SUMEL=SUMEL+IG*CI*BNEW*(1.-TEMP)*
     &		  BESSJ0(KT(IKT)*BNEW)
             IG=6-IG 
          ENDDO
          SUMEL=SUMEL*DBSM/3.
* add Coulomb amplitude
*----------------------
          TET=2.*ASIN(KT(IKT)/2./PCM)
          CALL FCOUL(PCM,ETA,TET,ARGAM,FC) 
          SUMEL=SUMEL+FC/PCM 	           
          SIGEL(IKT)=ABS(SUMEL)**2
        ENDDO
* loop over KY
*-------------
        DO IKY=0,NK
          PKY=KT(IKY)
* loop and integral over KX
*--------------------------
          SUMX=0.
          DO IKX=0,NK
            AUX=0.
            PKX=KT(IKX)
            PKT=SQRT(PKX*PKX+PKY*PKY)            
            IF(PKT.LT.KT(NK)) THEN
		     AUX=BETWEEN(PKT,SIGEL,KT,NK,NK)
              ELSE
               GOTO 89
            ENDIF 
89          SUMX=SUMX+SIMPK(IKX)*AUX
          ENDDO
          WRITE(20,100)PKY*HBARC,20.*SUMX/HBARC
        ENDDO
        PRINT*,' Output of elastic scattering X-section (mb/MeV/c) ',
     *	  FILENA
      ENDIF
*
90    FORMAT(1X,3(G11.5,2X)) 
*
100   FORMAT(1X,F10.4,2X,G11.5)
*
      RETURN
      END
*
*====================================================================
*
      SUBROUTINE COMPDIF(L,M,R,NR,B,NB,SNR,SNI,SCR,SCI,
     &    	WF,SIMPR,TERM1)
*
* Calculates a matrix used in the calculation of diffraction momentum 
* distributions.
*
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 R(0:NR),B(0:NB),SIMPR(0:NR),TERM1(0:NB),SNR(0:NB),
     &       IV(0:40),FAK(0:40),TERM2(0:NB),SCR(0:NB),SCI(0:NB),
     &       SNI(0:NB),FINT(0:NR),WF(0:NR)
      COMPLEX*16 CI,SUMPHI2,SUMBR2,SCT,SVT
      COMMON /COUL/ PCM,PCM2,ETA,ETA2
      COMMON /CONST1/ PI,FOUR_PI,HBARC,AMU,FINE,E2
      COMMON /RATIOS/ RATIOC,RATIOV
      COMMON /CONST2/ CI
      COMMON /GF1 / IV
      COMMON /GF2 / FAK
*
      NPHI=150
      DPHI=2.*PI/NPHI
* integral along z-axis
*----------------------
      CALL FINTEG(L,M,R,NR,WF,SIMPR,FINT)
* loop over B
*------------
      DO IB=0,NB
        BB=B(IB)	  
* integrate over R
*-----------------
        SUMBR0=0.
        SUMBR1=0.
        SUMBR2=0.
        DO IR=0,NR
          BR=R(IR)
          TOLD=TERM1(IB) 
* integrate over PHI
*-------------------
          SUMPHI0=0.
          SUMPHI1=0.
          SUMPHI2=0.
          SCT=0.
          DO IPHI=0,NPHI
            PHI=IPHI*DPHI
            CPHI=COS(PHI)
            BC=SQRT(ABS(BB**2+(RATIOC*BR)**2-2.*RATIOC*BB*BR*CPHI))
            BV=SQRT(ABS(BB**2+(RATIOV*BR)**2+2.*RATIOV*BB*BR*CPHI))
            IF(BC.LT.B(NB)) THEN
              SCT=BETWEEN(BC,SCR,B,NB,NB)+CI*BETWEEN(BC,SCI,B,NB,NB)
            ENDIF
            IF(BV.LT.B(NB)) THEN
              SVT=BETWEEN(BV,SNR,B,NB,NB)+CI*BETWEEN(BV,SNI,B,NB,NB)
            ENDIF
            SUMPHI1=SUMPHI1+FINT(IR)*ABS(1.-SCT*SVT)**2
            SUMPHI2=SUMPHI2+FINT(IR)*(1.-SCT*SVT)
            SUMPHI0=SUMPHI0+FINT(IR)
          ENDDO
*
          SUMBR0=SUMBR0+BR*SIMPR(IR)*SUMPHI0
          SUMBR1=SUMBR1+BR*SIMPR(IR)*SUMPHI1
          SUMBR2=SUMBR2+BR*SIMPR(IR)*SUMPHI2
        ENDDO
        TERM1(IB)=SUMBR1
        TERM2(IB)=ABS(SUMBR2)**2
* renormalize for accuracy
*-------------------------
        TERM1(IB)=TERM1(IB)/SUMBR0
        TERM2(IB)=TERM2(IB)/SUMBR0**2
* use completeness
*-----------------
        TERM1(IB)=TERM1(IB)-TERM2(IB)
      ENDDO
*
      RETURN
      END
*
*====================================================================
*
      SUBROUTINE GAUST(L,M,R,NR,KT,NK,BN,NB,NGS,SIMPR,WF,BLMSQ)
*
* Calculates the sum over the Gaussian expansion for d sigma / dp_y
*   #  p_y is one of the components of the transverse momentum  #
*
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NP=5,NGAUSS=40)
*
      REAL*8 R(0:NR),KT(0:NK),SIMPR(0:NR),BLMSQ(0:NK,0:NB),BN(0:NB),
     &     BETAN(0:NGAUSS),WF(0:NR),BETAN2(0:NGAUSS),BNN(0:NB),
     &     BNN2(0:NB),RHO(0:NR),RHO2(0:NR),Z(0:NR),Z2(0:NR),
     &     BMR2(0:NB,0:NR),TBR(0:NB,0:NR),BR2(0:NB,0:NR)
      COMPLEX*16	CI,SUM,SUMRHO,ALPHAC(NGAUSS),SG(0:NB,0:NR)
*
      COMMON /BT/ BETAN
      COMMON /GS2/ ALPHAC
      COMMON /COUL/ PCM,PCM2,ETA,ETA2
      COMMON /CONST2/ CI
*
      TOL=1.D-20
      DEN=0.39894
      DR=R(2)-R(1)
*
* repeated variables
*-------------------
      DO N=1,NGS
        BETAN2(N)=BETAN(N)*BETAN(N) 
      ENDDO
*
* mesh in transverse and longitudinal coordinates
*------------------------------------------------
      DO IR=0,NR
        RHO(IR)=R(IR)
        RHO2(IR)=RHO(IR)*RHO(IR) 
        Z(IR)=R(IR)	          
        Z2(IR)=Z(IR)*Z(IR) 
      ENDDO
*
* preliminary calculations
*-------------------------
      DO IB=0,NB
        BNN(IB)=BN(IB)
        BNN2(IB)=BNN(IB)*BNN(IB)
        DO IR=0,NR
          BMR2(IB,IR)=(BNN(IB)-RHO(IR))**2
          TBR(IB,IR)=2.*BNN(IB)*RHO(IR)
          BR2(IB,IR)=BNN2(IB)+RHO2(IR)
        ENDDO
      ENDDO
*
25    FORMAT(2X,'P=',I2)
26    FORMAT(4X,I2)
*
* sum over IP
*------------
      PRINT*,' Please be patient. This might take a while.'
      PRINT*,' Calculation ends after P=5   (P=-5, -4, ..., 4, 5)'
      PRINT*,' '
      PRINT*,' Starting:'
      SUMP=0.
      DO IP=-NP,NP
        IF(IP.EQ.-NP) THEN
            WRITE(*,25)IP	     
          ELSE
            WRITE(*,26)IP
        ENDIF
        MIP=IABS(M-IP) 
        IAP=IABS(IP)
* loops over BN and RHO
*---------------------- 
        DO IB=0,NB 
          DO IR=0,NR
* sum over gaussians
*-------------------
            SUM=0.
            IF(ETA.EQ.0.)SUM=BESSI(MIP,0.D0)
            DO N=1,NGS
              ARG=TBR(IB,IR)/BETAN2(N)
              IF(ARG.LT.350.) THEN
                 AUX=BESSI(MIP,ARG)
                 SUM=SUM+ALPHAC(N)*EXP(-BR2(IB,IR)/BETAN2(N))*AUX
                ELSE
                 MU=4.*MIP*MIP 
                 EA=8.*ARG
                 AUX=DEN/SQRT(ARG)*(1.-(MU-1.)/EA
     &     			   +(MU-1.)*(MU-9.)/2./EA/EA)
                 SUM=SUM+ALPHAC(N)*EXP(-BMR2(IB,IR)/BETAN2(N))*AUX
              ENDIF
            ENDDO
            SG(IB,IR)=SUM
* end sum over Gaussians
*----------------------- 
          ENDDO
        ENDDO
* 
* loop over KT
*-------------
        DO IKT=0,NK  
* loop over BN
*-------------
          DO IB=0,NB 
* integrate over Z
*-----------------
            SUMZ=0. 
            DO IZ=0,NR
* integrate over RHO
*-------------------
              SUMRHO=0.
              DO IR=0,NR
                WFMID=0.
                RMID=SQRT(Z2(IZ)+RHO2(IR))+TOL
                COST=Z(IZ)/RMID
                PLEG=PLGNDR(L,M,COST)
* interpolate function WF at point RMID
*--------------------------------------
                IF(RMID.LT.R(NR))WFMID=BETWEEN(RMID,WF,R,NR,NR)
* end interpolation
*------------------
                SUMRHO=SUMRHO+SIMPR(IR)*RHO(IR)*SG(IB,IR)*
     &   			BESSJ(IAP,KT(IKT)*RHO(IR))*WFMID*PLEG
              ENDDO
              SUMZ=SUMZ+SIMPR(IZ)*ABS(SUMRHO)**2          
            ENDDO
            BLMSQ(IKT,IB)=BLMSQ(IKT,IB)+2.*SUMZ
          ENDDO
        ENDDO
      ENDDO	           
* end integral over Z, RHO, and sum over P
*-----------------------------------------
*
      RETURN
      END
*
*====================================================================
*
      SUBROUTINE 
     &	GAUSLT(M,R,NR,KT,NK,BN,NB,SN,NGS,SIMPR,SIMPB,ZINT,DUM1)
*
* Calculates the sum over the Gaussian expansion
*
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NP=5,NGAUSS=40)
      REAL*8 R(0:NR),KT(0:NK),SIMPR(0:NR),DUM1(0:NK,0:NK),
     &  BN(0:NB),SN(0:NB),BETAN(0:NGAUSS),SIMPB(0:NB),BETAN2(0:NGAUSS),
     &     BNN(0:NB),BNN2(0:NB),RHO(0:NR),RHO2(0:NR),Z(0:NR),Z2(0:NR),
     &     BMR2(0:NB,0:NR),TBR(0:NB,0:NR),BR2(0:NB,0:NR)
      COMPLEX*16	CI,SUM,ZINT(0:NK,0:NR),SG(0:NB,0:NR),SUMRHO,
     &            ALPHAC(NGAUSS)
      COMMON /BT/ BETAN
      COMMON /GS2/ ALPHAC
      COMMON /COUL/ PCM,PCM2,ETA,ETA2
      COMMON /CONST2/ CI
*
* constants
*----------
      TOL=1.D-20
      DEN=0.39894  
      DR=R(2)-R(1)
*
* repeated variables
*-------------------
      DO N=1,NGS
        BETAN2(N)=BETAN(N)*BETAN(N) 
      ENDDO
*
* mesh in transverse and longitudinal coordinates
*------------------------------------------------
      DO IR=0,NR
        RHO(IR)=R(IR)
        RHO2(IR)=RHO(IR)*RHO(IR) 
        Z(IR)=R(IR)	          
        Z2(IR)=Z(IR)*Z(IR) 
      ENDDO
*
* preliminary calculations
*-------------------------
      DO IB=0,NB
        BNN(IB)=BN(IB)
        BNN2(IB)=BNN(IB)*BNN(IB)
        DO IR=0,NR
          BMR2(IB,IR)=(BNN(IB)-RHO(IR))**2
          TBR(IB,IR)=2.*BNN(IB)*RHO(IR)
          BR2(IB,IR)=BNN2(IB)+RHO2(IR)
        ENDDO
      ENDDO
*
25    FORMAT(2X,'P=',I2)
26    FORMAT(4X,I2)
*
* sum over IP
*------------
      PRINT*,' Please be patient. This might take a while.'
      PRINT*,' Calculation ends after P=5   (P=-5, -4, ..., 4, 5)'
      PRINT*,' '
      PRINT*,' Starting:'
      SUMP=0.
      DO IP=-NP,NP
        IF(IP.EQ.-NP) THEN
            WRITE(*,25)IP	     
          ELSE
            WRITE(*,26)IP
        ENDIF
        MIP=IABS(M-IP) 
        IAP=IABS(IP)
        DO IB=0,NB 
          DO IR=0,NR
* sum over Gaussians
*-------------------
            SUM=0.
            IF(ETA.EQ.0.)SUM=BESSI(MIP,0.D0)
            DO N=1,NGS
              ARG=TBR(IB,IR)/BETAN2(N)
              IF(ARG.LT.350.) THEN
                 AUX=BESSI(MIP,ARG)
                 SUM=SUM+ALPHAC(N)*
     &             EXP(-BR2(IB,IR)/BETAN2(N))*AUX
                ELSE
                 MU=4.*MIP*MIP 
                 EA=8.*ARG
                 AUX=DEN/SQRT(ARG)*(1.-(MU-1.)/EA
     &     			   +(MU-1.)*(MU-9.)/2./EA/EA)
                 SUM=SUM+ALPHAC(N)*EXP(-BMR2(IB,IR)/BETAN2(N))*AUX
              ENDIF
            ENDDO
* end sum over Gaussians 
*-----------------------
            SG(IB,IR)=SUM
          ENDDO
        ENDDO
* 
        DO IKT=0,NK  
          DO IKZ=0,NK
* integrate over BN     
*------------------   
            SUMB=0.  
            DO IB=0,NB 
* integrate over RHO
*-------------------
              SUMRHO=0.
              DO IR=0,NR
               SUMRHO=SUMRHO+SIMPR(IR)*RHO(IR)*SG(IB,IR)*
     &   			ZINT(IKZ,IR)*BESSJ(IAP,KT(IKT)*RHO(IR))
              ENDDO
              SUMB=SUMB+SIMPB(IB)*BN(IB)*(1.-SN(IB)**2)*ABS(SUMRHO)**2          
            ENDDO
            DUM1(IKZ,IKT)=DUM1(IKZ,IKT)+SUMB
          ENDDO
        ENDDO
      ENDDO	           
* end integral over Z, RHO, and sum over P
*-----------------------------------------
*
      RETURN
      END
*
*====================================================================
*
      SUBROUTINE FINTEG(L,M,R,NR,WF,SIMP,FINT)
*
* Calculates the integral over the z-coordinate of the function
*
*   WF(r) ** 2 * Plm[cos(theta)] 
*
* where WF(r) is a radial wavefunction, r = sqrt(b * b + z * z),
* cos(theta) = z / r, and Plm is the Legendre polynomial. 
*
* Output is a vector FINT(IR), where  B = IR * dR. 
*
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 R(0:NR),WF(0:NR),FINT(0:NR),SIMP(0:NR)
*
      TOL=1.D-20
      DR=R(2)-R(1)
      IAM=IABS(M)
*
      DO IR=0,NR
        B=R(IR)
* integrate over Z
*-----------------
        SUM=0.
        DO IZ=0,NR
          Z=IZ*DR
          RMID=SQRT(Z*Z+B*B)+TOL
          COST=Z/RMID
          PART1=PLGNDR(L,IAM,COST)
* interpolate function WF at point RMID
*--------------------------------------
          PART2=0.
          IF(RMID.LT.R(NR))PART2=BETWEEN(RMID,WF,R,NR,NR)
*
          SUM=SUM+2.*SIMP(IZ)*(PART1*PART2)**2
        ENDDO
* end integration over Z
*-----------------------
        FINT(IR)=SUM
      ENDDO	           
*
      RETURN
      END
*
*====================================================================
*
      SUBROUTINE ZINTEG(L,M,R,NR,WF,KZ,NK,SIMP,ZINT)
*
* Calculates the integral over the z-coordinate of the function
*
*   exp(-i * k_z * z) * WF(r) * Plm[cos(theta)]
*
* where WF(r) is a radial wavefunction, r = sqrt(b * b + z * z),
* cos(theta) = z / r, and Plm is the Legendre polynomial. 
*
* Output is a matrix ZINT(IKZ,IR), where K_z = IKZ * DK and B = IR * dR. 
*
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 R(0:NR),WF(0:NR),KZ(0:NK),SIMP(0:NR),IV(0:40)
      COMPLEX*16 CI,SUM,ZINT(0:NK,0:NR)
      COMMON /GF1 / IV
      COMMON /CONST2/ CI
*
      TOL=1.D-20
      DR=R(2)-R(1)
      DKZ=KZ(2)-KZ(1)
      IVLM=IV(L+M)
      IVMP=IVLM+1
      IVMM=IVLM-1
      IAM=IABS(M)
*
      DO IKZ=0,NK
        DO IR=0,NR
          B=R(IR)
* integrate over Z
*-----------------
          SUM=0.
          DO IZ=0,NR
            Z=IZ*DR
            PART0=COS(KZ(IKZ)*Z)
            PART1=SIN(KZ(IKZ)*Z)
            RMID=SQRT(Z*Z+B*B)+TOL
            COST=Z/RMID
            PART2=PLGNDR(L,IAM,COST)
* interpolate function WF at point RMID
*--------------------------------------
            PART3=0.
            IF(RMID.LT.R(NR))PART3=BETWEEN(RMID,WF,R,NR,NR)
*
            SUM=SUM+SIMP(IZ)*(IVMP*PART0+CI*IVMM*PART1)*PART2*PART3
          ENDDO
* end integration
*----------------
          ZINT(IKZ,IR)=SUM 
        ENDDO
      ENDDO	           
*
      RETURN
      END
*
*====================================================================
*
      SUBROUTINE PINTEG(R,NR,B,NB,S,SINT)
*
* Calculates the integral over phi of the function
*
*   | S ( sqrt(r^2 + b^2 - 2*r*b*cos(phi)) ) |^2
*
* where S is an S-matrix . 
*
* Output is a matrix PINT(IB,IR), where B = IB * DB and R = IR * dR. 
*
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NLG=300)
      REAL*8 R(0:NR),B(0:NB),SIMP(0:NLG),SINT(0:NB,0:NR),S(0:NB)
*
      COMMON /CONST1/ PI,FOUR_PI,HBARC,AMU,FINE,E2
*
      TOL=1.D-20
*
* Simpson coefficients for integral over PHI
*-------------------------------------------
      NPHI=150
      DPHI=2.*PI/NPHI
      SIMP(0)=1./3.*DPHI
      SIMP(NPHI)=1./3.*DPHI
      IS=4
      DO K=1,NPHI-1
        SIMP(K)=IS/3.*DPHI
        IS=6-IS
      ENDDO
*
      DR=R(2)-R(1)
      DB=B(2)-B(1)
*
      DO IB=0,NB
        DO IR=0,NR
* integrate over PHI
*-------------------
          SUM=0.
          DO IPHI=0,NPHI
             PHI=IPHI*DPHI
             BMID=SQRT(ABS(R(IR)**2+B(IB)**2-2.*R(IR)*B(IB)*COS(PHI)))
* interpolate function S
*-----------------------
             SFMID=0.
             IF(BMID.LT.B(NB))SFMID=BETWEEN(BMID,S,B,NB,NB)
*
             SUM=SUM+SIMP(IPHI)*SFMID**2
          ENDDO
* end integration
*----------------
          SINT(IB,IR)=SUM 
        ENDDO
      ENDDO	           
*
      RETURN
      END
*
*====================================================================
*
      SUBROUTINE COMPL(R,NR,NK,NB,SIMPR,ZINT,SINT,BLMSQ)
*
* Computes a matrix used in the calculation of longitudinal momentum 
* distributions.
*
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 R(0:NR),SIMPR(0:NR),BLMSQ(0:NK,0:NB),SINT(0:NB,0:NR) 
      COMPLEX*16	ZINT(0:NK,0:NR)
      COMMON /CONST1/ PI,FOUR_PI,HBARC,AMU,FINE,E2
*
* loops over BN and KZ 
*---------------------
      DO IB=0,NB 
        DO IKZ=0,NK
* integrate over RHO
*-------------------
          SUMRHO=0.
          DO IR=0,NR
            RHO=R(IR)
            SUMRHO=SUMRHO+SIMPR(IR)*RHO*ABS(ZINT(IKZ,IR))**2
     &           *SINT(IB,IR)
          ENDDO
          BLMSQ(IKZ,IB)=SUMRHO/2./PI
        ENDDO
      ENDDO	           
*
      RETURN
      END
*
*=====================================================================
*
      SUBROUTINE SMAT_FIT(IOPT,MA,BETA0,NDAT,AC)
*
* fits S-matrix (core+target) with an expansion in gaussians
*
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER*4 IA(MA)
      REAL*8 A1(MA),A2(MA),COVAR1(MA,MA),SIG1(NDAT),X(NDAT),Y1(NDAT),
     &	Y2(NDAT),AR(MA),AI(MA),SIG2(NDAT),COVAR2(MA,MA)
      COMPLEX*16 CI,AC(MA)
*
      EXTERNAL FUNCS
*
      COMMON/GF/AGF
      COMMON /COUL/ PCM,PCM2,ETA,ETA2
      COMMON /CONST2/ CI
*
      OPEN(30, FILE='SMAT_CORE_FIT.TXT', STATUS='UNKNOWN')
*
      AGF=BETA0
      NDATA=NDAT
*    
* read S-matrix (core+target) and add Coulomb phase
*--------------------------------------------------
      DO I=1,NDATA
        READ(22,*)X(I),Y1(I),Y2(I)
        CC=COS(2.*ETA*DLOG(PCM*X(I)))
        SS=SIN(2.*ETA*DLOG(PCM*X(I)))
        YR=Y1(I)
        YI=Y2(I)	   
        Y1(I)=YR*CC-YI*SS
        Y2(I)=YI*CC+YR*SS
        SIG1(I)=1
        SIG2(I)=1
        IF(ETA.EQ.0.)Y1(I)=Y1(I)-1. ! better fit if no Coulomb  
      ENDDO
      Y1O=Y1(1)
      Y2O=Y2(1)
*
* initialize input for fitting routines
*--------------------------------------
      IA(1)=0
      DO I=2,MA		        
	  IA(I)=1
      ENDDO	                 
*
* LFIT finds fitting coeficients and chi-squares
*-----------------------------------------------
      CALL LFIT(X,Y1,SIG1,NDAT,A1,IA,MA,COVAR1,MA,CHISQ1,FUNCS)
      CALL LFIT(X,Y2,SIG2,NDAT,A2,IA,MA,COVAR2,MA,CHISQ2,FUNCS)
*
      WRITE(30,*)'Chi-square deviation of the real and imaginary part'
      WRITE(30,*)' '
      WRITE(30,*)'CHISQR',CHISQ1
      WRITE(30,*)'CHISQI',CHISQ2
      WRITE(30,*)' '
*
*  These are the coefficients for the Gaussian expansion
*  F(X) = SUM_I A_I*EXP(-X**2/B_I**2), with  B_I=constant/I     
*----------------------------------------------------------
      WRITE(30,*)' Real and Imag. coefficients for the S-matrix fit'
      WRITE(30,*)' No. of coefficients = ',MA     
      WRITE(30,*)' '
      AR(1)=0.
      AI(1)=0.
      AC(1)=0.
      WRITE(10,*)AR(1),AI(1)
      DO I=2,MA
        AR(I)=A1(I)
        AI(I)=A2(I)
        WRITE(10,*)AR(I),AI(I)
        AC(I)=AR(I)+CI*AI(I)
      ENDDO
*
* output comparison between exact function and gaussian expansion
*----------------------------------------------------------------
      WRITE(30,*)' '
      WRITE(30,*)'Smat.  fit: 1st column = B (fm)'
      WRITE(30,*)'          : 2nd column = real part of S'
      WRITE(30,*)'          : 3rd column = Gaussian fit of real part'
      WRITE(30,*)'          : 4th column = imaginary part of S'
      WRITE(30,*)'          : 5th column = Gaussian fit of imag. part'
      WRITE(30,*)' '
      Y1(1)=Y1O
      Y2(1)=Y2O
      DO I=2,NDATA
        SUM1=0.
        SUM2=0.
        DO J=1,MA
          AGFJ=AGF/J
          SUM1=SUM1+AR(J)*EXP(-X(I)**2/AGFJ/AGFJ)
          SUM2=SUM2+AI(J)*EXP(-X(I)**2/AGFJ/AGFJ) 
        ENDDO 
        IF(ETA.NE.0.) THEN
            WRITE(30,100)X(I),Y1(I),SUM1,Y2(I),SUM2
           ELSE
            WRITE(30,100)X(I),Y1(I)+1.,SUM1+1.,Y2(I),SUM2
        ENDIF
      ENDDO
100   FORMAT(1X,5(G10.3,1X))
      PRINT*,' '
      IF(IOPT.NE.1.AND.IOPT.NE.5)
     &PRINT*,' Gauss fit of core+target S-matrix in SMAT_CORE_FIT.TXT'
      PRINT*,' '
*
      CLOSE(UNIT=22)
*		  	        			    
      END
*
*====================================================================
*
      SUBROUTINE FUNCS(X,AFUNC,MA)
*
* Gaussian basis with width decreasing with 1/N, and N_max = MA
*
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MMAX=40)
      REAL*8 AFUNC(MMAX)
      COMMON/GF/AGF
*
* basis functions are MA Gaussians, with widths 
* W = AGF/MA, AGF/(MA+1)...AGF
*-----------------------------
      DO I=1,MA
        AGFI=AGF/I
        AFUNC(I)=EXP(-X**2/AGFI/AGFI) 
      ENDDO
*
      RETURN
      END	 		        	           
*
*=====================================================================
*
      SUBROUTINE FCOUL(PCM,ETA,TET,ARGAM,FC)
*
*  elastic Coulomb amplitude
*
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 CI, FC
*
      COMMON /CONST2/ CI
*
      TOL=1.D-30
      FACT=-2.*ETA*DLOG(ABS(SIN((TET+TOL)/2.)))+2.*ARGAM
      SQRUT=ETA/(2.*PCM*SIN((TET+TOL)/2.)**2)
      FC=-EXP(CI*FACT)*SQRUT
*
      RETURN
      END
*
*=====================================================================
*
      SUBROUTINE CP0(ETA,ARGAM)
*
*  Coulomb phase
*
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (NEULER=5000)
      REAL*8 DELC(NEULER),DIF(NEULER)
*
      EULER=0.5772156649
*
* Coulomb phase
*--------------
	DELC(1)=ATAN(ETA)
	DO 49 L=2,NEULER
	  DELC(L)=DELC(L-1)+ATAN(ETA/FLOAT(L))
	  FASE=ETA*LOG(L+0.5)
	  DIF(L)=DELC(L)-FASE
	  DIFE=DIF(L)-DIF(L-1)
	  IF(ABS(DIFE).LT.1.0E-04)GO TO 490
   49 CONTINUE
  490 ARGAM=-DIF(L)
*
      RETURN
      END
*
*====================================================================
*
      SUBROUTINE SIGNNE(AP,ZP,AT,ZT,E0,IPAULI,SIGAV)
*
*  Isospin averaged nucleon-nucleon cross section (in mb) for lab 
*  energy E0
*
*  IPAULI=0: free nucleon-nucleon cross sections 
*  IPAULI=1: Pauli blocking corrected nucleon-nucleon cross sections
* 
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/PAUL/RHOEFF
*
      RHOEFF=0.17
      NP = NINT(AP-ZP)
      NT = NINT(AT-ZT)
*
* free nucleon nucleon cross section. Follows parametrization 
* of Charagi, Gupta, Phys. Rev. C 41, 1610 (1990)
*------------------------------------------------
      IF(IPAULI.EQ.0) THEN
* high energy fit
*----------------
        IF(E0.GT.10.D0) THEN
          BETA=SQRT(E0*E0 + 1863.D0*E0)/(E0+931.49386D0)
          SIGNP=-70.67D0-18.18D0/BETA+25.26D0/(BETA*BETA)+113.85D0*BETA
          SIGPP= 13.73D0-15.04D0/BETA+8.76D0/(BETA*BETA)+68.67D0
     & 		*(BETA**4)
        ELSE
*  low energy fit
*----------------
          BETA=SQRT(2.D0*E0/931.49386D0)
          SIGNP=((2.73/(0.35*E0+(1-0.0553*E0)**2))+
     &           (17.63/(6.8*E0+(1+0.344*E0)**2)))*1.0D3
          SIGPP=13.73-15.04/BETA+8.76/(BETA*BETA)+68.67*(BETA**4)
        ENDIF
      ENDIF
*
* including Pauli principle due to medium density RHOEFF. Uses 
* fit from Hussein, Rego, Bertulani, Phys. Rep. 201 (1991) 279
* and from Bertulani, J. Phys. G 27 (2001) L67
*---------------------------------------------
      IF(IPAULI.EQ.1) THEN
         BETA=SQRT(E0*E0 + 1863.D0*E0)/(E0+931.49386D0)
         DUM1=13.73-15.04/BETA+8.76/BETA**2+68.67*BETA**4
         DUM2=1.+7.772*E0**0.06*RHOEFF**1.48
         DUM3=1.+18.01*RHOEFF**1.46
         SIGPP=DUM1*DUM2/DUM3
         DUM1=-70.67-18.18/BETA+25.26/BETA**2+113.85*BETA
         DUM2=1.+20.88*E0**0.04*RHOEFF**2.02
         DUM3=1.+35.86*RHOEFF**1.9
         SIGNP=DUM1*DUM2/DUM3
      ENDIF
*
* isospin average
*----------------
      SIGAV=(NP*NT + ZP*ZT)*SIGPP + (NP*ZT + NT*ZP)*SIGNP
      SIGAV=SIGAV/(AP*AT)
*
* print np and pp cross sections if requested
*  PRINT*,'SIGNP =',SIGNP,' MB'
*  PRINT*,'SIGPP = SIGNN =',SIGPP,' MB'
*
      RETURN
      END
*
*====================================================================
*
      SUBROUTINE PHNNE(AP,ZP,AT,ZT,ECA,BNN,ALFANN)
*
* phase of nucleon-nucleon scattering amplitude from expanded table 1 
* of Hussein, Rego, Bertulani, Phys. Rep. 201 (1991) 279
* with aditional input from Lenzi, Vitturi, Zardi, Phys. Rev. C 40,
* 2114 (1989).
*
* The nucleon-nucleon scattering amplitude is given by
*     F(THETA ~ 0) = (K*SIGNN/4/PI)*(ALFANN+I)*EXP(-BNN*Q**2)
* where Q is the momentum transfer, K the nucleon momentum and SIGNN
* the nucleon-nucleon cross section.
* 
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (NDIM1=7,NDIM2=10)
	REAL*8 E1(NDIM1),E2(NDIM2),ALFP(NDIM2),ALFN(NDIM2),BP(NDIM2),
     &	BN(NDIM2),ALF(NDIM1)
	DATA E1/10.,30.,48.,40.,49.,85.,94./
      DATA ALF/0.8,0.87,0.89,0.9,0.94,1,1.07/
	DATA E2/100.,150.,200.,325.,425.,550.,650.,800.,1000.,2200./
	DATA ALFP/1.87,1.53,1.15,0.45,0.47,0.32,0.16,0.06,-0.09,
     &     	-0.17/
	DATA ALFN/1.,0.96,0.71,0.16,0.25,-0.24,-0.35,-0.2,
     &  	    -0.46,-0.5/
	DATA BP/0.66,0.57,0.56,0.26,0.21,0.04,0.07,0.09,0.09,0.12/
	DATA BN/0.36,0.58,0.68,0.36,0.27,0.085,0.09,0.12,0.12,0.14/
*
      NP = AP-ZP
      NT = AT-ZT
*
* interpolate
*------------
      IF(ECA.LT.10.OR.ECA.GT.2200.) THEN
	  PRINT*,'STOP: Energy/nucleon out of bounds, ECA=',eca
	  PRINT*,'      10. < ECA < 2200.'
        STOP
      ENDIF
*
      IF(ECA.GE.100.AND.ECA.LT.2200.)THEN
        ALFPP = BETWEEN(ECA,ALFP,E2,NDIM2,NDIM2) 
        ALFPN = BETWEEN(ECA,ALFN,E2,NDIM2,NDIM2)	   
        BPP   = BETWEEN(ECA,BP,E2,NDIM2,NDIM2) 
        BPN   = BETWEEN(ECA,BN,E2,NDIM2,NDIM2) 
*
* print n-n scattering amplitude parameters, if wanted
*      PRINT*,'ALFPP =',ALFPP,'ALFPN =',ALFPN
*      PRINT*,'BPP =',BPP,'BPN =',BPN
*      PRINT*
*
* isospin average
*----------------
        BNN=(NP*NT + ZP*ZT)*BPP + (NP*ZT + NT*ZP)*BPN
        BNN=BNN/(AP*AT)
        ALFANN=(NP*NT + ZP*ZT)*ALFPP + (NP*ZT + NT*ZP)*ALFPN
        ALFANN=ALFANN/(AP*AT)
	ENDIF
*
      IF(ECA.LT.100..AND.ECA.GT.10.) THEN
        ALFANN=BETWEEN(ECA,ALF,E1,NDIM1,NDIM1) 
        BNN=0.   
      ENDIF
*
      RETURN
      END
*
*======================================================================
*
      SUBROUTINE PHNUC(R,NR,WN,B,NB,PHASEN)
*
*  nuclear eikonal phase
*
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 B(NB),R(NR),Z(0:NR)
      COMPLEX*16 CI,PHASEN(NB),VNUCL,WN(NR)
      COMMON /CONST1/ PI,FOUR_PI,HBARC,AMU,FINE,E2
      COMMON/CONST2/CI
      COMMON/A1/ BETA,GAMMA
*
      DELR=R(3)-R(2)
      ZMAX=R(NR)
      ZMAX=ZMAX/NR
      DO I=0,NR
       Z(I)=I*DELR
      ENDDO
*
      DO K=1,NB
        PHASEN(K)=0.
*  integrate in z coordinate
*---------------------------
	  DO I=0,NR
          VNUC=0.
	    RV=SQRT(B(K)**2+Z(I)**2)
	    L=INT(RV/DELR)
	    AL=1.-(RV-DELR*L)/DELR
	    IF(L+1.LE.NR.AND.L.GE.1)VNUCL=AL*WN(L)+(1.-AL)*WN(L+1)
          PHASEN(K)=PHASEN(K)+VNUCL
	  ENDDO
	  PHASEN(K)=-2./HBARC/BETA*DELR*PHASEN(K)
      ENDDO
*
      RETURN
      END
*
*=====================================================================
*
      SUBROUTINE GAUSS(NRMAX,DR,DENS)
*
*  modified gaussian:  [1. + BETA * (R/ALPHA)**PWR] * EXP[-(R/ALPHA)**2]
*
*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL*8 DENS(NRMAX)
*
      WRITE(17,30)
   30 FORMAT(19H  Gaussian function)
*
      CALL SKIPCOM(12)      
      READ(12,*)ALPHA,BETA,PWR
      WRITE(17,60)ALPHA,BETA,PWR
   60 FORMAT(10H  ALPHA = ,F10.5,5X,7HBETA = ,F10.5,5X,6HPWR = ,
     &   F10.5,/)
*
      DO 100 I=1,NRMAX
        R=DR*DBLE(I)
  100   DENS(I)=(1.+BETA*(R/ALPHA)**PWR)*DEXP(-(R/ALPHA)**2)
*
      RETURN
      END
*
*=====================================================================
*
      SUBROUTINE YUKAWA(NRMAX,DR,DENS)
*
*  Yukawa:   R**PWR * EXP(-BETA*R/R) / BETA*R
*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL*8 DENS(NRMAX)
*
      WRITE(17,30)
   30 FORMAT(17H  Yukawa function)
      CALL SKIPCOM(12)      
      READ(12,*)BETA,PWR
      WRITE(17,60)BETA,PWR
   60 FORMAT(8HBETA =  ,F10.5,5X,7HPWR =   ,F10.5,/)
*
      DO 100 I=1,NRMAX
        R=DR*DBLE(I)
        BR=BETA*R
  100   DENS(I)=DEXP(-BR)/BR *R**PWR
*
      RETURN
      END
*
*=====================================================================
*
      SUBROUTINE FERMI(NRMAX,DR,DENS)
*
* 3 parameter fermi function
*   [1 + W * (R/AR)**PWR] / [1 + EXP((R-RR)/AR)]
*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL*8 DENS(NRMAX)
*
      WRITE(17,30)
   30 FORMAT(28H  3 Parameter Fermi function)
      CALL SKIPCOM(12)      
      READ(12,*)RR,AR,W,PWR
      WRITE(17,60)RR,AR,W,PWR
   60 FORMAT(5HRR = ,F10.5,5X,5HAR = ,F10.5,5X,4HW = ,F10.5,/)
*
      DO 100 I=1,NRMAX
        R=DR*DBLE(I)
        TEMP=(R-RR)/AR
        IF(TEMP.GT.30.)TEMP=30.
  100   DENS(I)=(1.0D0+W*(R/AR)**PWR)/(1.0D0+DEXP(TEMP))
*
      RETURN
      END
*
*=====================================================================
*
      SUBROUTINE SAXON(NRMAX,DR,DENS)
*
* Woods-saxon and derivative Woods-saxon functions raised to power
*  F = 1/[1 + EXP(R-RR)/AR)]**PWR   or  its  derivative   DF/DR 
*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL*8 DENS(NRMAX)
*
      CALL SKIPCOM(12)      
      READ(12,*)RR,AR,PWR,IDER
*
      IF(IDER.NE.0) GOTO 90
*
      WRITE(17,60)
   60 FORMAT(22H  Saxon-Woods function)
      WRITE(17,70)RR,AR,PWR
   70 FORMAT(5HRR = ,F10.5,5X,5HAR = ,F10.5,5X,6HPWR = ,F10.5,/)
*
      DO 80 I=1,NRMAX
        R=DR*DBLE(I)
        TEMP=(R-RR)/AR
        IF(TEMP.GT.30.)TEMP=30.
   80   DENS(I)=1./(1.0D0+DEXP(TEMP))**PWR
      RETURN
*
   90 WRITE(17,100)
  100 FORMAT(33H  Derivative Saxon-Woods function)
      WRITE(17,70)RR,AR,PWR
      PWRC=PWR/AR
      PWR1=PWR+1.0D0
      DO 110 I=1,NRMAX
        R=DR*DBLE(I)
        TEMP=(R-RR)/AR
        IF(TEMP.GT.30.)TEMP=30.
        EXPR=DEXP(TEMP)
  110   DENS(I)=-PWRC*EXPR/(1.D0+EXPR)**PWR1
*
      RETURN
      END
*
*=====================================================================
*
      SUBROUTINE DROP(AZ,AN,NRMAX,DR,DENS)
*
*  liquid drop density
*  ref.  W.D Myers, Nucl Phys A145 (1970) 387
*
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      REAL*8 DENS(NRMAX)
*
	PI=3.14159265359
      WRITE(17,30)
   30 FORMAT(21H  Liquid drop density,/)
      AKZ=0.75D0
      AKN=0.75D0
      IF(AZ.LT.0.0D0) AKZ=0.0D0
      IF(AN.LT.0.0D0) AKN=0.0D0
      AZ=DABS(AZ)
      AN=DABS(AN)
      AM=AZ+AN
      AI=(AN-AZ)/AM
*
      DELTA=(AI+0.0112D0*AZ**2/AM**1.66666666667D0)/
     *       (1.0D0+3.15D0/AM**0.333333333333D0)
      EPS=-0.147D0/AM**0.333333333333D0+0.33D0*DELTA**2
     *    +0.00248D0*AZ**2/AM**1.3333333333D0
      RRHO=1.16D0*AM**0.3333333333D0*(1.0D0+EPS)
      T=0.773D0*AM**0.3333333333D0*(AI-DELTA)
      RZ=RRHO-0.5D0*T
      RN=RRHO+0.5D0*T
      RZ=RZ-0.98D0/RZ
      RN=RN-0.98D0/RN
      A=0.55D0
      CZ=AKZ*AZ/(PI*RZ**3*(1.0D0+(PI*A/RZ)**2))
      CN=AKN*AN/(PI*RN**3*(1.0D0+(PI*A/RN)**2))
      VOLP=0.0D0
      VOLN=0.0D0
      RSQP=0.0D0
      RSQN=0.0D0
*
      DO 100 I=1,NRMAX
        R=DR*DBLE(I)
        RSQ=R*R
        R4=RSQ*RSQ
        TEMP=(R-RZ)/A
        IF(TEMP.GT.30.)TEMP=30.
        TP=CZ/(1.0D0+DEXP(TEMP))
        TEMP=(R-RN)/A
        IF(TEMP.GT.30.)TEMP=30.
        TN=CN/(1.0D0+DEXP(TEMP))
        VOLP=VOLP+TP*RSQ
        VOLN=VOLN+TN*RSQ
        RSQP=RSQP+R4*TP
        RSQN=RSQN+R4*TN
        DENS(I)=TP+TN
  100 CONTINUE
*
   80 RETURN
      END
*
*====================================================================
*
      SUBROUTINE LFIT(X,Y,SIG,NDAT,A,IA,MA,COVAR,NPC,CHISQ,FUNCS)
* 
* Given a set of data points X[1..NDAT], Y[1..NDAT] with individual
* standard deviations SIG[1..NDAT], use chi-square minimization to fit
* for some or all of the coefficients A[1..MA] of a function that depends
* linearly on A:
*     y = Sum_i a_i x afunc_i(x).
* The input array IA[1..MA] indicates by nonzero entries those components of
* A that should be fitted for, and by zero entries those components that should
* be held fixed at their input values. The program returns values for A[1..MA],
* chi-square = CHISQ, and the covariance matrix COVAR[1..MA]. (Parameters held
* fixed will return zero covariances.) The user supplies a routine FUNCS(X,AFUNC,MA)
* that return the MA basis functions evaluated at x = X in the array AFUNC[1..MA].
*
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER*4 MA,IA(MA),NPC,NDAT,MMAX
      REAL*8 CHISQ,A(MA),COVAR(NPC,NPC),SIG(NDAT),X(NDAT),Y(NDAT)
      EXTERNAL FUNCS
      PARAMETER (MMAX=50)
*
* Uses COVSRT, GAUSSJ
*--------------------
      INTEGER*4 I,J,K,L,M,MFIT
      REAL*8 SIG2I,SUM,WT,YM,AFUNC(MMAX),BETA(MMAX)
*
      MFIT=0
      DO 11 J=1,MA
        IF(IA(J).NE.0) MFIT=MFIT+1
11    CONTINUE
*
      IF(MFIT.EQ.0) PAUSE 'LFIT: NO PARAMETERS TO BE FITTED'
      DO 13 J=1,MFIT
        DO 12 K=1,MFIT
          COVAR(J,K)=0.
12      CONTINUE
        BETA(J)=0.
13    CONTINUE
*
      DO 17 I=1,NDAT
        CALL FUNCS(X(I),AFUNC,MA)
        YM=Y(I)
        IF(MFIT.LT.MA) THEN
          DO 14 J=1,MA
            IF(IA(J).EQ.0) YM=YM-A(J)*AFUNC(J)
14        CONTINUE
        ENDIF
        SIG2I=1./SIG(I)**2
        J=0
        DO 16 L=1,MA
          IF (IA(L).NE.0) THEN
            J=J+1
            WT=AFUNC(L)*SIG2I
            K=0
            DO 15 M=1,L
              IF (IA(M).NE.0) THEN
                K=K+1
                COVAR(J,K)=COVAR(J,K)+WT*AFUNC(M)
              ENDIF
15          CONTINUE
            BETA(J)=BETA(J)+YM*WT
          ENDIF
16      CONTINUE
17    CONTINUE
*
      DO 19 J=2,MFIT
        DO 18 K=1,J-1
          COVAR(K,J)=COVAR(J,K)
18      CONTINUE
19    CONTINUE
*
      CALL GAUSSJ(COVAR,MFIT,NPC,BETA,1,1)
      J=0
      DO 21 L=1,MA
        IF(IA(L).NE.0) THEN
          J=J+1
          A(L)=BETA(J)
        ENDIF
21    CONTINUE
      CHISQ=0.
*
      DO 23 I=1,NDAT
        CALL FUNCS(X(I),AFUNC,MA)
        SUM=0.
        DO 22 J=1,MA
          SUM=SUM+A(J)*AFUNC(J)
22      CONTINUE
        CHISQ=CHISQ+((Y(I)-SUM)/SIG(I))**2
23    CONTINUE
      CALL COVSRT(COVAR,NPC,MA,IA,MFIT)
*
      RETURN
      END
*
*====================================================================
*
      SUBROUTINE COVSRT(COVAR,NPC,MA,IA,MFIT)
*
* Expand in storage the covariance matrix COVAR, so as to take into 
* account parameters that are being held fixed. (For the latter, return
* zero covariances).
*
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER*4 MA,MFIT,NPC,IA(MA)
      REAL*8 COVAR(NPC,NPC)
      INTEGER*4 I,J,K
      REAL*8 SWAP
*
      DO 12 I=MFIT+1,MA
        DO 11 J=1,I
          COVAR(I,J)=0.
          COVAR(J,I)=0.
11      CONTINUE
12    CONTINUE
      K=MFIT
*
      DO 15 J=MA,1,-1
        IF(IA(J).NE.0)THEN
          DO 13 I=1,MA
            SWAP=COVAR(I,K)
            COVAR(I,K)=COVAR(I,J)
            COVAR(I,J)=SWAP
13        CONTINUE
          DO 14 I=1,MA
            SWAP=COVAR(K,I)
            COVAR(K,I)=COVAR(J,I)
            COVAR(J,I)=SWAP
14        CONTINUE
          K=K-1
        ENDIF
15    CONTINUE
*
      RETURN
      END
*
*====================================================================
*
      SUBROUTINE GAUSSJ(A,N,NP,B,M,MP)
*
* Linear equation solution by Gauss-Jordan elimination (se Eq. 2.1.1
* of Numerical Recipes). A[1..N][1..N] is the input matrix. B[1..N][1..M]
* is the input containing the M right-hand side vectors. On output, A is 
* replaced by its matrix inverse, and B is replaced by the corresponding 
* set of solution vectors. 
*
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER*4 M,MP,N,NP,NMAX
      REAL*8 A(NP,NP),B(NP,MP)
      PARAMETER (NMAX=50)
      INTEGER*4 I,ICOL,IROW,J,K,L,LL,INDXC(NMAX),INDXR(NMAX),IPIV(NMAX)
      REAL*8 BIG,DUM,PIVINV
*
      DO 11 J=1,N
        IPIV(J)=0
11    CONTINUE
*
      DO 22 I=1,N
        BIG=0.
        DO 13 J=1,N
          IF(IPIV(J).NE.1)THEN
            DO 12 K=1,N
              IF (IPIV(K).EQ.0) THEN
                IF (ABS(A(J,K)).GE.BIG)THEN
                  BIG=ABS(A(J,K))
                  IROW=J
                  ICOL=K
                ENDIF
              ELSE IF (IPIV(K).GT.1) THEN
                PAUSE 'SINGULAR MATRIX IN GAUSSJ'
              ENDIF
12          CONTINUE
          ENDIF
13      CONTINUE
        IPIV(ICOL)=IPIV(ICOL)+1
        IF (IROW.NE.ICOL) THEN
          DO 14 L=1,N
            DUM=A(IROW,L)
            A(IROW,L)=A(ICOL,L)
            A(ICOL,L)=DUM
14        CONTINUE
          DO 15 L=1,M
            DUM=B(IROW,L)
            B(IROW,L)=B(ICOL,L)
            B(ICOL,L)=DUM
15        CONTINUE
        ENDIF
        INDXR(I)=IROW
        INDXC(I)=ICOL
        IF (A(ICOL,ICOL).EQ.0.) PAUSE 'SINGULAR MATRIX IN GAUSSJ'
        PIVINV=1./A(ICOL,ICOL)
        A(ICOL,ICOL)=1.
        DO 16 L=1,N
          A(ICOL,L)=A(ICOL,L)*PIVINV
16      CONTINUE
        DO 17 L=1,M
          B(ICOL,L)=B(ICOL,L)*PIVINV
17      CONTINUE
        DO 21 LL=1,N
          IF(LL.NE.ICOL)THEN
            DUM=A(LL,ICOL)
            A(LL,ICOL)=0.
            DO 18 L=1,N
              A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
18          CONTINUE
            DO 19 L=1,M
              B(LL,L)=B(LL,L)-B(ICOL,L)*DUM
19          CONTINUE
          ENDIF
21      CONTINUE
22    CONTINUE
*
      DO 24 L=N,1,-1
        IF(INDXR(L).NE.INDXC(L))THEN
          DO 23 K=1,N
            DUM=A(K,INDXR(L))
            A(K,INDXR(L))=A(K,INDXC(L))
            A(K,INDXC(L))=DUM
23        CONTINUE
        ENDIF
24    CONTINUE
*
      RETURN
      END
*
*====================================================================
*
      REAL*8 FUNCTION PLGNDR(L,M,X)
*
* Calculates the Legendre polynomial of order L,M. For X < 0, use
*      PLM(-X) = (-1)**(L+M) * PLM(X).
*
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER*4 L,M
      REAL*8 X
      INTEGER*4 I,LL
      REAL*8 FACT,PLL,PMM,PMMP1,SOMX2
*
      IF(M.LT.0.OR.M.GT.L.OR.ABS(X).GT.1.)PAUSE
     *'BAD ARGUMENTS IN PLGNDR'
*
      PMM=1.
      IF(M.GT.0) THEN
        SOMX2=SQRT((1.-X)*(1.+X))
        FACT=1.
        DO 11 I=1,M
          PMM=-PMM*FACT*SOMX2
          FACT=FACT+2.
11      CONTINUE
      ENDIF
*
      IF(L.EQ.M) THEN
        PLGNDR=PMM
      ELSE
        PMMP1=X*(2*M+1)*PMM
        IF(L.EQ.M+1) THEN
          PLGNDR=PMMP1
        ELSE
          DO 12 LL=M+2,L
            PLL=(X*(2*LL-1)*PMMP1-(LL+M-1)*PMM)/(LL-M)
            PMM=PMMP1
            PMMP1=PLL
12        CONTINUE
          PLGNDR=PLL
        ENDIF
      ENDIF
*
      RETURN
      END
*
*=====================================================================
*      
      REAL*8 FUNCTION BESSJ(N,X)
*
* Bessel function J of order N. Do not use for N < 0. Remember that 
*    J(-N,X) = (-)**N * J(N,X).
*
      IMPLICIT REAL*8(A-H,O-Z)
      INTEGER*4 N,IACC
      REAL*8 X,BIGNO,BIGNI
      PARAMETER (IACC=40,BIGNO=1.E10,BIGNI=1.E-10)
*
* uses BESSJ0,BESSJ1
*-------------------
      INTEGER*4 J,JSUM,M
      REAL*8 AX,BJ,BJM,BJP,SUM,TOX,BESSJ0,BESSJ1
*
      IF(N.EQ.0) THEN
         BESSJ=BESSJ0(X)
         RETURN
      ENDIF	   	
      IF(N.EQ.1) THEN
         BESSJ=BESSJ1(X)
         IF(X.LT.0.)BESSJ=-BESSJ
         RETURN
      ENDIF	   	
*
      AX=ABS(X)
      IF(AX.EQ.0.)THEN
        BESSJ=0.
      ELSE IF(AX.GT.FLOAT(N))THEN
        TOX=2./AX
        BJM=BESSJ0(AX)
        BJ=BESSJ1(AX)
        DO 11 J=1,N-1
          BJP=J*TOX*BJ-BJM
          BJM=BJ
          BJ=BJP
11      CONTINUE
        BESSJ=BJ
      ELSE
        TOX=2./AX
        M=2*((N+INT(SQRT(FLOAT(IACC*N))))/2)
        BESSJ=0.
        JSUM=0
        SUM=0.
        BJP=0.
        BJ=1.
        DO 12 J=M,1,-1
          BJM=J*TOX*BJ-BJP
          BJP=BJ
          BJ=BJM
          IF(ABS(BJ).GT.BIGNO)THEN
            BJ=BJ*BIGNI
            BJP=BJP*BIGNI
            BESSJ=BESSJ*BIGNI
            SUM=SUM*BIGNI
          ENDIF
          IF(JSUM.NE.0)SUM=SUM+BJ
          JSUM=1-JSUM
          IF(J.EQ.N)BESSJ=BJP
12      CONTINUE
        SUM=2.*SUM-BJ
        BESSJ=BESSJ/SUM
      ENDIF
      IF(X.LT.0..AND.MOD(N,2).EQ.1)BESSJ=-BESSJ
*
      RETURN
      END
*
*=====================================================================
*      
      REAL*8 FUNCTION BESSJ0(X)
*
*  Bessel function of order 0
*
      IMPLICIT REAL*8(A-H,O-Z)
      DATA P1,P2,P3,P4,P5/1.D0,-.1098628627D-2,.2734510407D-4,
     * -.2073370639D-5,.2093887211D-6/,Q1,Q2,Q3,Q4,Q5/-.1562499995D-1,
     * .1430488765D-3,-.6911147651D-5,.7621095161D-6,-.934945152D-7/
      DATA R1,R2,R3,R4,R5,R6/57568490574.D0,-13362590354.D0,
     * 651619640.7D0,-11214424.18D0,77392.33017D0,-184.9052456D0/,
     * S1,S2,S3,S4,S5,S6/57568490411.D0,1029532985.D0,
     * 9494680.718D0,59272.64853D0,267.8532712D0,1.D0/
*
      IF(ABS(X).LT.8.)THEN
      Y=X**2
      BESSJ0=(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))
     * /(S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6)))))
      ELSE
*
      AX=ABS(X)
      Z=8./AX
      Y=Z**2
      XX=AX-.785398164
      BESSJ0=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+
     * Y*P5))))-Z*SIN(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)))))
      ENDIF
*
      RETURN
      END
*
*=====================================================================
*      
      REAL*8 FUNCTION BESSJ1(X)
*
*  Bessel function of order 1
*
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 X
      REAL*8 AX,XX,Z
      REAL*8  P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6,
     *S1,S2,S3,S4,S5,S6,Y
      SAVE P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6,S1,S2,S3,S4,
     *S5,S6
      DATA R1,R2,R3,R4,R5,R6/72362614232.D0,-7895059235.D0,
     *242396853.1D0,-2972611.439D0,15704.48260D0,-30.16036606D0/,S1,S2,
     *S3,S4,S5,S6/144725228442.D0,2300535178.D0,18583304.74D0,
     *99447.43394D0,376.9991397D0,1.D0/
      DATA P1,P2,P3,P4,P5/1.D0,.183105D-2,-.3516396496D-4,
     *.2457520174D-5,-.240337019D-6/, Q1,Q2,Q3,Q4,Q5/.04687499995D0,
     *-.2002690873D-3,.8449199096D-5,-.88228987D-6,.105787412D-6/
*
      ONE=1.
      IF(ABS(X).LT.8.)THEN
        Y=X**2
        BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))/(S1+Y*(S2+Y*(S3+
     *Y*(S4+Y*(S5+Y*S6)))))
      ELSE
        AX=ABS(X)
        Z=8./AX
        Y=Z**2
        XX=AX-2.356194491
        BESSJ1=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y*
     *P5))))-Z*SIN(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)))))*SIGN(ONE,X)
      ENDIF
*
      RETURN
      END
*
*=====================================================================
*      
      REAL*8 FUNCTION BESSI(N,X)
*
* Modified Bessel function I of order N.  Do not use for N < 0. Remember 
* that I(-N,X) = I(N,X).
*
      IMPLICIT REAL*8(A-H,O-Z)
      INTEGER*4 N,IACC
      REAL*8 X,BIGNO,BIGNI
      PARAMETER (IACC=40,BIGNO=1.0E10,BIGNI=1.0E-10)
* uses BESSI0
      INTEGER*4 J,M
      REAL*8 BI,BIM,BIP,TOX,BESSI0
*
      IF(N.EQ.0) THEN
         BESSI=BESSI0(X)
         RETURN
      ENDIF	   	
      IF(N.EQ.1) THEN
         BESSI=BESSI1(X)
         IF (X.LT.0.) BESSI=-BESSI
         RETURN
      ENDIF	   	
*
      IF (X.EQ.0.) THEN
        BESSI=0.
      ELSE
        TOX=2.0/ABS(X)
        BIP=0.0
        BI=1.0
        BESSI=0.
        M=2*((N+INT(SQRT(FLOAT(IACC*N)))))
        DO 11 J=M,1,-1
          BIM=BIP+FLOAT(J)*TOX*BI
          BIP=BI
          BI=BIM
          IF (ABS(BI).GT.BIGNO) THEN
            BESSI=BESSI*BIGNI
            BI=BI*BIGNI
            BIP=BIP*BIGNI
          ENDIF
          IF (J.EQ.N) BESSI=BIP
11      CONTINUE
        BESSI=BESSI*BESSI0(X)/BI
        IF (X.LT.0..AND.MOD(N,2).EQ.1) BESSI=-BESSI
      ENDIF
*
      RETURN
      END
*
*=====================================================================
*      
      REAL*8 FUNCTION BESSI0(X)
*
* Modified Bessel function I of order 0.
*
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 AX
      REAL*8 P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,Y
      SAVE P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9
      DATA P1,P2,P3,P4,P5,P6,P7/1.0D0,3.5156229D0,3.0899424D0,
     *1.2067492D0,0.2659732D0,0.360768D-1,0.45813D-2/
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,0.1328592D-1,
     *0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1,0.2635537D-1,
     *-0.1647633D-1,0.392377D-2/
*
      IF (ABS(X).LT.3.75) THEN
        Y=(X/3.75)**2
        BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
      ELSE
        AX=ABS(X)
        Y=3.75/AX
        BESSI0=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*
     *(Q7+Y*(Q8+Y*Q9))))))))
      ENDIF
*
      RETURN
      END
*
*=====================================================================
*      
      REAL*8 FUNCTION BESSI1(X)
*
* Modified Bessel function I of order 1.
*
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 X
      REAL*8 AX
      REAL*8 P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,Y
      SAVE P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9
      DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0,
     *0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,-0.3988024D-1,
     *-0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1,-0.2895312D-1,
     *0.1787654D-1,-0.420059D-2/
*
      IF (ABS(X).LT.3.75) THEN
        Y=(X/3.75)**2
        BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
      ELSE
        AX=ABS(X)
        Y=3.75/AX
        BESSI1=(EXP(AX)/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*
     *(Q7+Y*(Q8+Y*Q9))))))))
        IF(X.LT.0.)BESSI1=-BESSI1
      ENDIF
*
      RETURN
      END
*
*=====================================================================
*      
      SUBROUTINE GFV(IE)
*
*  Calculates sign and factorials  of integers and half int.
*
*     IV(N)  =  (-1)**N
*     FAK(N) =  N!
*     FAD(N) =  N!!
*
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 IV(0:40),FAK(0:40),FAD(0:40)
      COMMON /GF1 / IV
      COMMON /GF2 / FAK
      COMMON /GF3 / FAD
*
      IF (IE.GT.40) THEN
       WRITE(12,*) '******    STOP IN GFV: IE LARGER 40 '
       STOP
      ENDIF
*
      IV(0)  = +1
      FAK(0) =  1.D0
      FAD(0)=1.
      FAD(1)=1.
      DO 10 I=1,IE
         IV(I)  = -IV(I-1)
10       FAK(I) = I*FAK(I-1)
      DO 11 I=3,IE,2
         FAD(I) = I*FAD(I-2)
11    CONTINUE 
      DO 12 I=2,IE,2
         FAD(I) = I*FAD(I-2)
12    CONTINUE 
*
      RETURN
	END
*
*====================================================================
*
      REAL*8 FUNCTION BETWEEN(R,FUNCT,RGRID,NPTS,NDIM)
*
*  Calculates, by interpolation, the value of a real function at an 
*  arbitrary point R, when the value of the function (stored in array FUNCT) 
*  is known on a grid of points. The values of the NPTS points at which the 
*  function is known is stored in array RGRID.
*  NDIM is the externally defined dimensions of the arrays FUNCT and RGRID
*
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 FUNCT(0:NDIM),RGRID(0:NDIM)
*
      DO 30 K=1,NPTS
      NST=0
      IF(RGRID(K).LT.R) GOTO 30
      NST=MAX0(K-3,1)
      GOTO 33
   30 CONTINUE
*
   33 IF((NST.GT.NPTS-5).OR.(NST.EQ.0)) NST=NPTS-5
      X1=RGRID(NST+0)
      X2=RGRID(NST+1)
      X3=RGRID(NST+2)
      X4=RGRID(NST+3)
      X5=RGRID(NST+4)
      X6=RGRID(NST+5)
      Y1=FUNCT(NST+0)
      Y2=FUNCT(NST+1)
      Y3=FUNCT(NST+2)
      Y4=FUNCT(NST+3)
      Y5=FUNCT(NST+4)
      Y6=FUNCT(NST+5)
      PII1=(X1-X2)*(X1-X3)*(X1-X4)*(X1-X5)*(X1-X6)
      PII2=(X2-X1)*(X2-X3)*(X2-X4)*(X2-X5)*(X2-X6)
      PII3=(X3-X1)*(X3-X2)*(X3-X4)*(X3-X5)*(X3-X6)
      PII4=(X4-X1)*(X4-X2)*(X4-X3)*(X4-X5)*(X4-X6)
      PII5=(X5-X1)*(X5-X2)*(X5-X3)*(X5-X4)*(X5-X6)
      PII6=(X6-X1)*(X6-X2)*(X6-X3)*(X6-X4)*(X6-X5)
  777 XD1=R-X1
      XD2=R-X2
      XD3=R-X3
      XD4=R-X4
      XD5=R-X5
      XD6=R-X6
      PI1= XD2*XD3*XD4*XD5*XD6
      PI2= XD1*XD3*XD4*XD5*XD6
      PI3= XD1*XD2*XD4*XD5*XD6
      PI4= XD1*XD2*XD3*XD5*XD6
      PI5= XD1*XD2*XD3*XD4*XD6
      PI6= XD1*XD2*XD3*XD4*XD5
      BETWEEN=Y1*PI1/PII1+Y2*PI2/PII2+Y3*PI3/PII3+Y4*PI4/PII4+
     + Y5*PI5/PII5+Y6*PI6/PII6
*
      RETURN
      END
*
*=====================================================================
*
      SUBROUTINE FINFOUT(FILENA,NRMAX,DR,FOUT)
*
*  Reads function FIN, defined on a X-mesh with NIN equally-spaced 
*  points and returns function FOUT, defined on a mesh with step size 
*  DR and NRMAX points.
* The first row in FIN must be the number of rows (NIN). The following 
* rows are the R x FIN input data.
* 
* FIN   = function read in from cards
* FOUT  = interpolated function returned by subroutine
*
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(NBIG=10000)
      REAL*8 FOUT(NRMAX),FIN(0:NBIG),RIN(0:NBIG)
      CHARACTER*12 FILENA
      OPEN(UNIT=13,FILE=FILENA,STATUS='UNKNOWN')
*
      READ(13,*)NIN
      READ(13,*)(RIN(I),FIN(I),I=1,NIN)
      DRIN=RIN(3)-RIN(2)
*
      IF(DR.EQ.DRIN.AND.NIN.EQ.NRMAX) GOTO 165
*
      DO J=1,NRMAX
        RX=DR*J
        FOUT(J)=BETWEEN(RX,FIN,RIN,NIN,NBIG)
      ENDDO
      GO TO 200
*
  165 DO J=1,NIN
        FOUT(J)=FIN(J)
      ENDDO
*
  200 RETURN
      END
*
*=====================================================================
*
      SUBROUTINE SKIPCOM(IUNIT)
*
* allows comments in input files
*
      CHARACTER*1 ICOM
*
1     READ(IUNIT,'(A1)')ICOM
      IF(ICOM.EQ.'*')GOTO 1
      BACKSPACE IUNIT
*
      RETURN
      END
*
*=====================================================================
*
      SUBROUTINE RSIMP(FA,SUM,H,NN)
*
*  Simpson's integral of discrete function  fa (real)
*  sum=result (real), h=integration  step
*  nn= maximum number of points (< 10000)
*
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NBIG=10000)
      REAL*8 FA(NN),S,SUM                        
*
      L=NN                                     
      S=FA(1)-FA(L)                           
       DO 1 I=2,L,2                            
 1     S=S+4.D0*FA(I)+2.D0*FA(I+1)           
       SUM=H*S/3.D0                       
*
       RETURN                                  
       END
*
*=====================================================================
*
      SUBROUTINE FOURIER0(FUNC,R,NR,FOUR,Q,NQ)
*
*  Computes the Fourier transform of a spherically symmetric function
*  FUNC, defined on a grid R with NR points. Returns Fourier
*  transform in function FOUR on a grid Q with NQ points
*
*  based on Filon's formula in M Abramowitz and I A Stegun (1965)
*  Handbook of Mathematical Functions, Dover, New York.
*
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 K,KH,KH2,KH3,FUNC(NR),FOUR(NQ),R(NR),Q(NQ)
*
      DR=R(3)-R(2)
      DQ=Q(3)-Q(2)
      PI=3.14159265
      PIPI=2.0D0*PI
      CONST=4.0D0*PI
*
      DO 100 J=1,NQ
        K=DQ*DBLE(J)
        KH=K*DR
        KH2=KH*KH
        KH3=KH2*KH
        SINKH=DSIN(DMOD(KH,PIPI))
        COSKH=DCOS(DMOD(KH,PIPI))
        A=1.0D0/KH+SINKH*COSKH/KH2-2.0D0*SINKH*SINKH/KH3
        B=(1.0D0+COSKH*COSKH)/KH2-2.0D0*SINKH*COSKH/KH3
        C=4.0D0*(SINKH/KH3-COSKH/KH2)
        EVEN=0.0D0
        ODD=0.0
        DO 20 I=2,NR,2
          F=DBLE(I)
          X=F*KH
          RR=F*DR
   20   EVEN=EVEN+FUNC(I)*RR*DSIN(DMOD(X,PIPI))
        DO 25 I=1,NR,2
          F=DBLE(I)
          X=F*KH
          RR=F*DR
   25   ODD=ODD+FUNC(I)*RR*DSIN(DMOD(X,PIPI))
        RR=DR*DBLE(NR)
        X=RR*K
        AA=FUNC(NR)*RR
        S=-AA*DSIN(DMOD(X,PIPI))+2.0D0*EVEN
  100 FOUR(J)=CONST*DR*(-A*AA*DCOS(DMOD(X,PIPI))+B*S+C*ODD)/K
*
      RETURN
      END
*
*======================================================================
*
      SUBROUTINE NUCNAME(IZ,IA,ANAME)
*
* returns the chemical symbol as a 5 character string in the form nnnaa
*
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER ELEMENT(112)*2,NN(11)*1,ANAME*5
      DATA ELEMENT/
     *'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne','Na','Mg','Al',
     *'Si','P ','S ','Cl','Ar','K ','Ca','Sc','Ti','V ','Cr','Mn','Fe',
     *'Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr','Rb','Sr','Y ',
     *'Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te',
     *'I ','Xe','Cs','Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb',
     *'Dy','Ho','Er','Tm','Yb','Lu','Hf','Ta','Wr','Re','Os','Ir','Pt',
     *'Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra','Ac','Th','Pa',
     *'U ','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No','Lr','Rf',
     *'Db','Sg','Bh','Hs','Mt','10','11','12'/
      DATA NN/'0','1','2','3','4','5','6','7','8','9',' '/
*
	IF(IZ.GT.0.AND.IZ.LT.112)THEN
        DO 10 I=1,2
          NA=IA
          IA=IA/10
          NA=NA-IA*10
          IF(NA.EQ.0.AND.IA.EQ.0) NA=10
          ANAME(4-I:4-I)=NN(NA+1)
 10     CONTINUE
*
        IF(IA.EQ.0) IA=10
        ANAME(1:1)=NN(IA+1)
        ANAME(4:5)=ELEMENT(IZ)
*
	ELSE IF(IZ.EQ.0.)THEN
	  ANAME='n'
	ELSE
	  ANAME='   '
	END IF
*
      RETURN
      END
*
*======================================================================
*
      FUNCTION EXCESS(ZZ,AA,APX)
*
*  Nuclear excess mass table 
*
*  Excess returns the mass excess (in MeV) from Audi & Wapstra, if possible. 
*  If the nucleus is not in Audi & Wapstra, it tries to extrapolate from the 
*  table. If this too fails, the routine uses a liquid drop fit to the  table. 
*  Note that the liquid drop fit is always used when the liquid drop pairing 
*  parameter apx is less than or equal to zero.
*
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NUMZ=112,NUMISO=2931)
*
      DIMENSION Z(3,3),ZM(3),C(9)
*
      COMMON/XCESS/IZ(NUMZ),IA(NUMZ),P(NUMISO)
      DATA C/-2.308,-7.263,30.496,19.344,52.151,0.311,0.720,
     1        2.034,11.236/
*
* if APX < 0 use liquid-drop
*---------------------------
      APR=APX
      IF(APX.LT.0.) APR=C(9)
      IF(APX.LE.0.) GO TO 110
*
* neutron
*--------
      IF(ZZ.EQ.0) THEN
        EXCESS=P(1)
        RETURN
      ENDIF
*
* nuclei
*-------
      IZETA=ZZ+0.5
      IAA=AA+0.5
      EXCESS=0.
      IERR=0
*
      IND1=IAA-IA(IZETA)
      INDEX=IZ(IZETA)+IND1
      IF(INDEX.GT.NUMISO) GO TO 110
      IF(INDEX.GE.IZ(IZETA+1).OR.IND1.LT.0) GO TO 20
      EXCESS=P(INDEX)
*
      RETURN
*
* Bad luck. Interpolate mass excess
*----------------------------------
 20   IZED=IZETA
      IZD=2
      IF(IND1.LT.0) IZD=-2
      Z(1,1)=0.
      Z(2,1)=0.
      Z(3,1)=0.
      Z(3,2)=0.
      Z(3,3)=0.
      ZM(1)=0.
      ZM(2)=0.
      ZM(3)=0.
 30   IZED=IZED+IZD
      IF(IZED.LT.8.OR.IZED.GE.NUMZ) GO TO 50
      IND1=IAA-IA(IZED)
      INDEX=IZ(IZED)+IND1
      IF(INDEX.LT.IZ(IZED+1).AND.IND1.GE.0) GO TO 40
      IF(Z(1,1).GT.0.OR.INDEX.GT.NUMISO) GO TO 50
      GO TO 30
 40   Z(1,1)=Z(1,1)+1.
      Z(2,1)=Z(2,1)+IZED
      Z(3,1)=Z(3,1)+IZED**2
      Z(3,2)=Z(3,2)+IZED**3
      Z(3,3)=Z(3,3)+IZED**4
      EX=P(INDEX)
      ZM(1)=ZM(1)+EX
      ZM(2)=ZM(2)+EX*IZED
      ZM(3)=ZM(3)+EX*IZED**2
      GO TO 30
 50   IF(Z(1,1).LT.3) GO TO 110
      Z(1,2)=Z(2,1)
      Z(2,2)=Z(3,1)
      Z(1,3)=Z(3,1)
      Z(2,3)=Z(3,2)
      DO 80 I=1,2
      IF(ABS(Z(I,I)).LT.1.0E-35) GO TO 110
      I1=I+1
      DO 70 J=I1,3
      FAC=Z(J,I)/Z(I,I)
      ZM(J)=ZM(J)-FAC*ZM(I)
      DO 60 K=I1,3
      Z(J,K)=Z(J,K)-FAC*Z(I,K)
 60   CONTINUE
 70   CONTINUE
 80   CONTINUE
      IF(ABS(Z(3,3)).LT.1.0E-35) GO TO 110
      ZM(3)=ZM(3)/Z(3,3)
      EXCESS=ZM(3)
      DO 100 I=2,1,-1
      I1=I+1
      DO 90 J=I1,3
      ZM(I)=ZM(I)-Z(I,J)*ZM(J)
 90   CONTINUE
      ZM(I)=ZM(I)/Z(I,I)
      EXCESS=ZM(I)+IZETA*EXCESS
 100  CONTINUE
      IERR=1
*
      RETURN
*
* liquid-drop mass excess
*------------------------
 110  A3=AA**(-1./3.)
      D=(1.-2.*ZZ/AA)/(1.+C(8)*A3)
      D=D*D
      EXCESS=ZZ*C(1)+AA*(C(2)+C(3)*D+A3*(C(4)+C(5)*D+C(6)*A3))
      EXCESS=EXCESS+C(7)*ZZ*(ZZ-1.)*A3
      D=2.*(1+INT(ZZ+.5)/2+INT(AA-ZZ+.5)/2)-AA
      EXCESS=EXCESS-APR*D/SQRT(AA)
      IERR=-1
*
      RETURN
      END
*
*======================================================================
*
      BLOCK DATA
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NUMZ=112,NUMISO=2931)
*
* Wapstra and Audi mass table
      COMMON/XCESS/IZ(NUMZ),IA(NUMZ),P(NUMISO)
*
      DATA P(   1),P(   2),P(   3)/    8.07132,    7.28897,   13.13572/
      DATA P(   4),P(   5),P(   6)/   14.94979,   25.92778,   36.83398/
      DATA P(   7),P(   8),P(   9)/   41.86376,   14.93120,    2.42491/
      DATA P(  10),P(  11),P(  12)/   11.38624,   17.59412,   26.11027/
      DATA P(  13),P(  14),P(  15)/   31.59798,   40.81837,   48.81001/
      DATA P(  16),P(  17),P(  18)/   25.32018,   11.67888,   14.08631/
      DATA P(  19),P(  20),P(  21)/   14.90767,   20.94620,   24.95390/
      DATA P(  22),P(  23),P(  24)/   33.05023,   40.79586,   50.09600/
      DATA P(  25),P(  26),P(  27)/   37.99600,   18.37447,   15.76949/
      DATA P(  28),P(  29),P(  30)/    4.94166,   11.34758,   12.60658/
      DATA P(  31),P(  32),P(  33)/   20.17397,   25.07640,   33.65845/
      DATA P(  34),P(  35),P(  36)/   39.88240,   27.86786,   22.92100/
      DATA P(  37),P(  38),P(  39)/   12.41570,   12.05076,    8.66798/
      DATA P(  40),P(  41),P(  42)/   13.36890,   16.56221,   23.66373/
      DATA P(  43),P(  44),P(  45)/   28.96694,   37.08169,   43.71631/
      DATA P(  46),P(  47),P(  48)/   52.32200,   59.36400,   35.09406/
      DATA P(  49),P(  50),P(  51)/   28.91365,   15.69857,   10.65053/
      DATA P(  52),P(  53),P(  54)/     .00000,    3.12501,    3.01989/
      DATA P(  55),P(  56),P(  57)/    9.87314,   13.69412,   21.03659/
      DATA P(  58),P(  59),P(  60)/   24.92404,   32.83339,   37.56006/
      DATA P(  61),P(  62),P(  63)/   45.96000,   52.58300,   39.69900/
      DATA P(  64),P(  65),P(  66)/   24.96052,   17.33808,    5.34546/
      DATA P(  67),P(  68),P(  69)/    2.86342,     .10144,    5.68343/
      DATA P(  70),P(  71),P(  72)/    7.87082,   13.11714,   15.86045/
      DATA P(  73),P(  74),P(  75)/   21.76649,   25.23191,   32.08089/
      DATA P(  76),P(  77),P(  78)/   37.73500,   47.04000,   32.04784/
      DATA P(  79),P(  80),P(  81)/   23.11073,    8.00646,    2.85539/
      DATA P(  82),P(  83),P(  84)/   -4.73700,    -.80900,    -.78206/
      DATA P(  85),P(  86),P(  87)/    3.33356,    3.79691,    8.06174/
      DATA P(  88),P(  89),P(  90)/    9.28435,   14.61637,   18.97446/
      DATA P(  91),P(  92),P(  93)/   27.14400,   35.16400,   33.60800/
      DATA P(  94),P(  95),P(  96)/   16.77700,   10.68026,    1.95170/
      DATA P(  97),P(  98),P(  99)/     .87343,   -1.48741,    -.01740/
      DATA P( 100),P( 101),P( 102)/    -.04758,    2.79378,    3.32952/
      DATA P( 103),P( 104),P( 105)/    7.54451,   11.26638,   18.28816/
      DATA P( 106),P( 107),P( 108)/   25.05003,   33.22600,   40.29600/
      DATA P( 109),P( 110),P( 111)/   23.99240,   16.48517,    5.30678/
      DATA P( 112),P( 113),P( 114)/    1.75106,   -7.04193,   -5.73172/
      DATA P( 115),P( 116),P( 117)/   -8.02434,   -5.15364,   -5.94752/
      DATA P( 118),P( 119),P( 120)/   -2.05870,     .42988,    7.09351/
      DATA P( 121),P( 122),P( 123)/   11.27859,   18.02059,   22.23662/
      DATA P( 124),P( 125),P( 126)/   30.84200,   37.17600,   25.31800/
      DATA P( 127),P( 128),P( 129)/   12.92862,    6.84486,   -2.18426/
      DATA P( 130),P( 131),P( 132)/   -5.18210,   -9.52949,   -8.41760/
      DATA P( 133),P( 134),P( 135)/   -9.35746,   -6.90247,   -5.58086/
      DATA P( 136),P( 137),P( 138)/   -1.03358,    2.61871,    8.59442/
      DATA P( 139),P( 140),P( 141)/   12.66376,   18.30366,   25.50989/
      DATA P( 142),P( 143),P( 144)/   32.50900,   41.15300,   17.57053/
      DATA P( 145),P( 146),P( 147)/   10.91168,    -.39677,   -5.47267/
      DATA P( 148),P( 149),P( 150)/  -13.93338,  -13.19273,  -16.21448/
      DATA P( 151),P( 152),P( 153)/  -14.58650,  -15.01875,  -10.66119/
      DATA P( 154),P( 155),P( 156)/   -8.88223,   -3.21509,    -.79560/
      DATA P( 157),P( 158),P( 159)/    5.20423,    8.45092,   16.29200/
      DATA P( 160),P( 161),P( 162)/   20.91200,   29.10000,   26.11900/
      DATA P( 163),P( 164),P( 165)/   18.18300,    6.76721,    -.05504/
      DATA P( 166),P( 167),P( 168)/   -8.91574,  -12.21034,  -17.19683/
      DATA P( 169),P( 170),P( 171)/  -16.85055,  -18.21551,  -15.87237/
      DATA P( 172),P( 173),P( 174)/  -14.95418,  -11.06207,   -8.50492/
      DATA P( 175),P( 176),P( 177)/   -2.86224,    -.05808,    5.91638/
      DATA P( 178),P( 179),P( 180)/    9.60370,   15.74200,   20.40000/
      DATA P( 181),P( 182),P( 183)/   32.16400,   23.77200,   10.75476/
      DATA P( 184),P( 185),P( 186)/    3.82531,   -7.14462,  -12.38443/
      DATA P( 187),P( 188),P( 189)/  -21.49279,  -21.89503,  -24.43288/
      DATA P( 190),P( 191),P( 192)/  -22.94896,  -24.08086,  -20.49238/
      DATA P( 193),P( 194),P( 195)/  -19.95656,  -14.35976,  -12.40064/
      DATA P( 196),P( 197),P( 198)/   -6.52419,   -3.74461,    2.14200/
      DATA P( 199),P( 200),P( 201)/    5.40300,   11.83000,   14.99700/
      DATA P( 202),P( 203),P( 204)/   31.99700,   18.87200,   10.97300/
      DATA P( 205),P( 206),P( 207)/    -.75298,   -7.16102,  -16.95191/
      DATA P( 208),P( 209),P( 210)/  -20.20056,  -24.44099,  -24.30532/
      DATA P( 211),P( 212),P( 213)/  -26.33773,  -24.55755,  -24.85761/
      DATA P( 214),P( 215),P( 216)/  -20.25084,  -18.99471,  -14.46610/
      DATA P( 217),P( 218),P( 219)/  -12.64969,   -8.33687,   -4.84377/
      DATA P( 220),P( 221),P( 222)/     .08400,    3.08300,    9.20300/
      DATA P( 223),P( 224),P( 225)/   14.10300,   22.19700,   25.97000/
      DATA P( 226),P( 227),P( 228)/   17.50700,    4.07311,   -3.15888/
      DATA P( 229),P( 230),P( 231)/  -14.06281,  -19.04493,  -26.01598/
      DATA P( 232),P( 233),P( 234)/  -26.58624,  -29.93185,  -28.84637/
      DATA P( 235),P( 236),P( 237)/  -30.66396,  -26.89622,  -26.86108/
      DATA P( 238),P( 239),P( 240)/  -23.16134,  -22.84955,  -18.60193/
      DATA P( 241),P( 242),P( 243)/  -17.24195,  -12.48202,  -10.88000/
      DATA P( 244),P( 245),P( 246)/   -4.82500,    -.40100,    7.09800/
      DATA P( 247),P( 248),P( 249)/   12.10000,   20.50200,   26.55700/
      DATA P( 250),P( 251),P( 252)/   13.14300,    4.44300,   -7.06444/
      DATA P( 253),P( 254),P( 255)/  -13.33070,  -21.00351,  -24.44057/
      DATA P( 256),P( 257),P( 258)/  -29.01351,  -29.52189,  -31.76152/
      DATA P( 259),P( 260),P( 261)/  -29.79798,  -29.80065,  -27.55773/
      DATA P( 262),P( 263),P( 264)/  -27.33915,  -24.98733,  -24.02939/
      DATA P( 265),P( 266),P( 267)/  -19.99106,  -18.90932,  -14.79200/
      DATA P( 268),P( 269),P( 270)/  -11.22500,   -4.79700,    -.10200/
      DATA P( 271),P( 272),P( 273)/    7.20000,   12.60300,   20.08300/
      DATA P( 274),P( 275),P( 276)/   11.29600,   -2.17908,   -9.38134/
      DATA P( 277),P( 278),P( 279)/  -18.37827,  -23.04821,  -30.23044/
      DATA P( 280),P( 281),P( 282)/  -30.94803,  -34.71476,  -33.24184/
      DATA P( 283),P( 284),P( 285)/  -35.03989,  -33.06726,  -34.42207/
      DATA P( 286),P( 287),P( 288)/  -31.97753,  -32.26204,  -29.71933/
      DATA P( 289),P( 290),P( 291)/  -29.72074,  -25.90835,  -23.22200/
      DATA P( 292),P( 293),P( 294)/  -16.59900,  -13.09700,   -6.29700/
      DATA P( 295),P( 296),P( 297)/   -1.70500,    5.80000,   20.41800/
      DATA P( 298),P( 299),P( 300)/    6.76300,   -1.48100,  -11.16711/
      DATA P( 301),P( 302),P( 303)/  -17.42508,  -24.79924,  -28.80169/
      DATA P( 304),P( 305),P( 306)/  -33.80684,  -33.53502,  -35.55887/
      DATA P( 307),P( 308),P( 309)/  -35.02132,  -36.59304,  -35.81021/
      DATA P( 310),P( 311),P( 312)/  -36.60803,  -35.41893,  -35.69689/
      DATA P( 313),P( 314),P( 315)/  -32.12448,  -30.32005,  -25.35263/
      DATA P( 316),P( 317),P( 318)/  -22.00200,  -16.19900,  -11.99800/
      DATA P( 319),P( 320),P( 321)/   -5.59800,    -.57000,   13.15300/
      DATA P( 322),P( 323),P( 324)/    4.43900,   -6.43920,  -13.16061/
      DATA P( 325),P( 326),P( 327)/  -22.05905,  -27.27626,  -34.84611/
      DATA P( 328),P( 329),P( 330)/  -35.13749,  -38.54677,  -38.40844/
      DATA P( 331),P( 332),P( 333)/  -41.46909,  -40.81253,  -43.13491/
      DATA P( 334),P( 335),P( 336)/  -42.33970,  -44.21474,  -41.29005/
      DATA P( 337),P( 338),P( 339)/  -39.57146,  -35.88651,  -32.50914/
      DATA P( 340),P( 341),P( 342)/  -27.89800,  -23.58500,  -18.11800/
      DATA P( 343),P( 344),P( 345)/  -13.23700,   -7.12000,   13.89800/
      DATA P( 346),P( 347),P( 348)/    2.84100,   -4.93700,  -14.16801/
      DATA P( 349),P( 350),P( 351)/  -20.52639,  -28.64221,  -32.12093/
      DATA P( 352),P( 353),P( 354)/  -36.18762,  -37.81581,  -41.06934/
      DATA P( 355),P( 356),P( 357)/  -41.75864,  -44.33163,  -44.49281/
      DATA P( 358),P( 359),P( 360)/  -46.55228,  -44.53751,  -43.21880/
      DATA P( 361),P( 362),P( 363)/  -40.38026,  -37.96800,  -34.46527/
      DATA P( 364),P( 365),P( 366)/  -30.33900,  -25.46700,  -21.38700/
      DATA P( 367),P( 368),P( 369)/  -15.77000,  -11.14000,    9.10100/
      DATA P( 370),P( 371),P( 372)/    1.23200,   -8.85021,  -15.71300/
      DATA P( 373),P( 374),P( 375)/  -25.12088,  -29.32032,  -37.54830/
      DATA P( 376),P( 377),P( 378)/  -39.00691,  -44.12534,  -44.93173/
      DATA P( 379),P( 380),P( 381)/  -48.48700,  -48.55804,  -51.42585/
      DATA P( 382),P( 383),P( 384)/  -49.72685,  -49.46403,  -46.82461/
      DATA P( 385),P( 386),P( 387)/  -45.76429,  -41.80545,  -39.13206/
      DATA P( 388),P( 389),P( 390)/  -34.55800,  -31.56800,  -26.11900/
      DATA P( 391),P( 392),P( 393)/  -22.69100,  -16.75000,   10.33000/
      DATA P( 394),P( 395),P( 396)/    -.24200,   -8.16900,  -18.02400/
      DATA P( 397),P( 398),P( 399)/  -23.84600,  -31.87359,  -37.07393/
      DATA P( 400),P( 401),P( 402)/  -42.00393,  -44.47466,  -47.95618/
      DATA P( 403),P( 404),P( 405)/  -49.21754,  -52.19749,  -51.43741/
      DATA P( 406),P( 407),P( 408)/  -51.84461,  -49.88673,  -49.14730/
      DATA P( 409),P( 410),P( 411)/  -46.23936,  -44.37637,  -40.38026/
      DATA P( 412),P( 413),P( 414)/  -37.91180,  -33.06803,  -30.35700/
      DATA P( 415),P( 416),P( 417)/  -25.02000,  -21.65700,    5.99000/
      DATA P( 418),P( 419),P( 420)/   -2.13600,  -13.53500,  -19.41200/
      DATA P( 421),P( 422),P( 423)/  -29.47093,  -34.55236,  -42.81532/
      DATA P( 424),P( 425),P( 426)/  -45.32544,  -50.25446,  -51.44476/
      DATA P( 427),P( 428),P( 429)/  -55.41280,  -55.28064,  -56.92832/
      DATA P( 430),P( 431),P( 432)/  -55.10330,  -55.28860,  -52.39299/
      DATA P( 433),P( 434),P( 435)/  -51.93078,  -47.85084,  -46.82620/
      DATA P( 436),P( 437),P( 438)/  -42.76488,  -41.17203,  -35.52700/
      DATA P( 439),P( 440),P( 441)/  -33.34700,  -27.60000,    6.39900/
      DATA P( 442),P( 443),P( 444)/   -5.11400,  -12.37000,  -22.26300/
      DATA P( 445),P( 446),P( 447)/  -28.99700,  -37.61054,  -42.62147/
      DATA P( 448),P( 449),P( 450)/  -48.23696,  -50.70114,  -54.68363/
      DATA P( 451),P( 452),P( 453)/  -55.55127,  -57.70639,  -56.90555/
      DATA P( 454),P( 455),P( 456)/  -57.48486,  -55.90226,  -55.47310/
      DATA P( 457),P( 458),P( 459)/  -52.91444,  -51.73517,  -48.46563/
      DATA P( 460),P( 461),P( 462)/  -46.75168,  -43.10022,  -40.89258/
      DATA P( 463),P( 464),P( 465)/  -36.49600,  -33.70100,   13.56300/
      DATA P( 466),P( 467),P( 468)/     .75500,   -6.62300,  -18.10800/
      DATA P( 469),P( 470),P( 471)/  -24.58200,  -34.47150,  -40.21731/
      DATA P( 472),P( 473),P( 474)/  -48.32914,  -50.94128,  -56.24841/
      DATA P( 475),P( 476),P( 477)/  -57.47501,  -60.60101,  -60.17571/
      DATA P( 478),P( 479),P( 480)/  -62.14885,  -60.65842,  -61.40693/
      DATA P( 481),P( 482),P( 483)/  -58.91749,  -58.89790,  -55.77931/
      DATA P( 484),P( 485),P( 486)/  -55.07923,  -51.28805,  -50.31930/
      DATA P( 487),P( 488),P( 489)/  -46.57469,  -44.23700,  -39.40200/
      DATA P( 490),P( 491),P( 492)/    1.63900,   -9.57600,  -17.19500/
      DATA P( 493),P( 494),P( 495)/  -27.27400,  -33.91600,  -42.63914/
      DATA P( 496),P( 497),P( 498)/  -48.00533,  -54.02371,  -56.03501/
      DATA P( 499),P( 500),P( 501)/  -59.33967,  -59.84143,  -62.22361/
      DATA P( 502),P( 503),P( 504)/  -61.64422,  -62.89505,  -61.42810/
      DATA P( 505),P( 506),P( 507)/  -61.83702,  -59.78931,  -59.16423/
      DATA P( 508),P( 509),P( 510)/  -56.05226,  -55.32143,  -51.82832/
      DATA P( 511),P( 512),P( 513)/  -51.04586,  -46.75200,  -44.96300/
      DATA P( 514),P( 515),P( 516)/  -40.60400,   -3.79100,  -11.43900/
      DATA P( 517),P( 518),P( 519)/  -22.65400,  -29.37900,  -39.20610/
      DATA P( 520),P( 521),P( 522)/  -45.32991,  -53.89965,  -56.07548/
      DATA P( 523),P( 524),P( 525)/  -60.22302,  -61.15113,  -64.46810/
      DATA P( 526),P( 527),P( 528)/  -64.21677,  -66.74269,  -65.50922/
      DATA P( 529),P( 530),P( 531)/  -67.09590,  -65.12259,  -66.02873/
      DATA P( 532),P( 533),P( 534)/  -63.74246,  -63.48603,  -60.37772/
      DATA P( 535),P( 536),P( 537)/  -59.48520,  -55.88964,  -54.67870/
      DATA P( 538),P( 539),P( 540)/  -50.22600,  -48.52200,  -43.80800/
      DATA P( 541),P( 542),P( 543)/  -41.61000,  -36.48700,  -33.72000/
      DATA P( 544),P( 545),P( 546)/   -2.62700,  -13.46000,  -21.69400/
      DATA P( 547),P( 548),P( 549)/  -31.62400,  -38.60100,  -47.30527/
      DATA P( 550),P( 551),P( 552)/  -51.65997,  -56.35155,  -58.34121/
      DATA P( 553),P( 554),P( 555)/  -61.97957,  -62.79452,  -65.57616/
      DATA P( 556),P( 557),P( 558)/  -65.42081,  -67.25972,  -66.25433/
      DATA P( 559),P( 560),P( 561)/  -67.30016,  -65.54189,  -65.73991/
      DATA P( 562),P( 563),P( 564)/  -62.96034,  -62.76423,  -60.06300/
      DATA P( 565),P( 566),P( 567)/  -59.15900,  -55.70300,  -54.30600/
      DATA P( 568),P( 569),P( 570)/  -50.31000,  -48.48400,  -43.95700/
      DATA P( 571),P( 572),P( 573)/  -41.65600,  -35.49900,   -6.56700/
      DATA P( 574),P( 575),P( 576)/  -14.92300,  -25.72800,  -32.68600/
      DATA P( 577),P( 578),P( 579)/  -42.29311,  -47.25741,  -54.18311/
      DATA P( 580),P( 581),P( 582)/  -56.34243,  -61.16735,  -62.20930/
      DATA P( 583),P( 584),P( 585)/  -65.99953,  -65.90778,  -68.89631/
      DATA P( 586),P( 587),P( 588)/  -67.87716,  -70.00404,  -68.41493/
      DATA P( 589),P( 590),P( 591)/  -69.55943,  -67.32168,  -68.12842/
      DATA P( 592),P( 593),P( 594)/  -65.41000,  -65.70920,  -62.46842/
      DATA P( 595),P( 596),P( 597)/  -62.04289,  -58.60414,  -57.22206/
      DATA P( 598),P( 599),P( 600)/  -53.39800,  -51.77735,  -46.12800/
      DATA P( 601),P( 602),P( 603)/  -42.06600,   -4.74100,  -15.90100/
      DATA P( 604),P( 605),P( 606)/  -23.98600,  -34.12100,  -39.99800/
      DATA P( 607),P( 608),P( 609)/  -47.34800,  -51.99635,  -56.68930/
      DATA P( 610),P( 611),P( 612)/  -58.83473,  -62.65291,  -63.72130/
      DATA P( 613),P( 614),P( 615)/  -66.87669,  -67.08293,  -69.32092/
      DATA P( 616),P( 617),P( 618)/  -68.90471,  -70.13683,  -68.58650/
      DATA P( 619),P( 620),P( 621)/  -69.70385,  -68.05401,  -68.46420/
      DATA P( 622),P( 623),P( 624)/  -66.20290,  -65.87415,  -63.66207/
      DATA P( 625),P( 626),P( 627)/  -62.48799,  -59.06775,  -57.98275/
      DATA P( 628),P( 629),P( 630)/  -52.94600,  -49.49000,  -44.39500/
      DATA P( 631),P( 632),P( 633)/   -8.37400,  -17.00000,  -27.76800/
      DATA P( 634),P( 635),P( 636)/  -33.72900,  -42.24300,  -46.91000/
      DATA P( 637),P( 638),P( 639)/  -54.42473,  -56.41056,  -61.62130/
      DATA P( 640),P( 641),P( 642)/  -62.65376,  -66.97696,  -67.09364/
      DATA P( 643),P( 644),P( 645)/  -70.56033,  -69.90490,  -72.58556/
      DATA P( 646),P( 647),P( 648)/  -71.29713,  -73.42201,  -71.85591/
      DATA P( 649),P( 650),P( 651)/  -73.21289,  -71.21414,  -71.86207/
      DATA P( 652),P( 653),P( 654)/  -69.48800,  -69.44775,  -66.30274/
      DATA P( 655),P( 656),P( 657)/  -65.62344,  -61.00400,  -58.39500/
      DATA P( 658),P( 659),P( 660)/  -53.38400,  -50.04900,   -6.39900/
      DATA P( 661),P( 662),P( 663)/  -18.05200,  -24.96400,  -33.82300/
      DATA P( 664),P( 665),P( 666)/  -39.52100,  -47.05600,  -51.82100/
      DATA P( 667),P( 668),P( 669)/  -56.64376,  -58.87696,  -63.08062/
      DATA P( 670),P( 671),P( 672)/  -64.34032,  -67.89219,  -68.22946/
      DATA P( 673),P( 674),P( 675)/  -70.95628,  -70.85960,  -73.03246/
      DATA P( 676),P( 677),P( 678)/  -72.28958,  -73.91618,  -72.81621/
      DATA P( 679),P( 680),P( 681)/  -73.63599,  -72.11797,  -72.53275/
      DATA P( 682),P( 683),P( 684)/  -70.32344,  -69.88009,  -66.08000/
      DATA P( 685),P( 686),P( 687)/  -63.51900,  -59.40100,  -56.28100/
      DATA P( 688),P( 689),P( 690)/  -51.64200,  -47.29200,  -32.91900/
      DATA P( 691),P( 692),P( 693)/  -41.72200,  -46.49100,  -54.14800/
      DATA P( 694),P( 695),P( 696)/  -56.29748,  -61.94000,  -63.09200/
      DATA P( 697),P( 698),P( 699)/  -67.89443,  -68.21628,  -72.21262/
      DATA P( 700),P( 701),P( 702)/  -72.16882,  -75.25156,  -74.59905/
      DATA P( 703),P( 704),P( 705)/  -77.02567,  -75.91694,  -77.75941/
      DATA P( 706),P( 707),P( 708)/  -76.38908,  -77.59344,  -75.34009/
      DATA P( 709),P( 710),P( 711)/  -75.94980,  -72.42860,  -70.54095/
      DATA P( 712),P( 713),P( 714)/  -66.58249,  -63.87814,  -59.59700/
      DATA P( 715),P( 716),P( 717)/  -56.43000,  -50.88800,  -47.19900/
      DATA P( 718),P( 719),P( 720)/  -32.79800,  -38.89200,  -46.40900/
      DATA P( 721),P( 722),P( 723)/  -51.59000,  -56.59200,  -59.15277/
      DATA P( 724),P( 725),P( 726)/  -63.53264,  -65.30596,  -69.13882/
      DATA P( 727),P( 728),P( 729)/  -70.28869,  -73.23393,  -73.45190/
      DATA P( 730),P( 731),P( 732)/  -76.06798,  -75.88885,  -77.97437/
      DATA P( 733),P( 734),P( 735)/  -77.49595,  -79.00911,  -77.77631/
      DATA P( 736),P( 737),P( 738)/  -78.61060,  -75.63995,  -73.85749/
      DATA P( 739),P( 740),P( 741)/  -70.73215,  -68.56982,  -64.61308/
      DATA P( 742),P( 743),P( 744)/  -61.51092,  -56.58336,  -53.00200/
      DATA P( 745),P( 746),P( 747)/  -47.80400,  -32.30400,  -40.97600/
      DATA P( 748),P( 749),P( 750)/  -46.10000,  -54.11277,  -56.88529/
      DATA P( 751),P( 752),P( 753)/  -62.16956,  -64.24160,  -68.97871/
      DATA P( 754),P( 755),P( 756)/  -70.17141,  -74.15971,  -74.44221/
      DATA P( 757),P( 758),P( 759)/  -77.89335,  -77.69365,  -80.58857/
      DATA P( 760),P( 761),P( 762)/  -79.98184,  -82.43104,  -81.48061/
      DATA P( 763),P( 764),P( 765)/  -83.26595,  -80.70999,  -79.69215/
      DATA P( 766),P( 767),P( 768)/  -76.72482,  -74.96309,  -71.31293/
      DATA P( 769),P( 770),P( 771)/  -68.78826,  -64.02600,  -61.14100/
      DATA P( 772),P( 773),P( 774)/  -56.03900,  -53.03000,  -47.91600/
      DATA P( 775),P( 776),P( 777)/  -32.30400,  -38.11700,  -46.23400/
      DATA P( 778),P( 779),P( 780)/  -51.72551,  -57.22242,  -60.48055/
      DATA P( 781),P( 782),P( 783)/  -64.82583,  -66.93577,  -70.79659/
      DATA P( 784),P( 785),P( 786)/  -72.17278,  -75.45644,  -76.18903/
      DATA P( 787),P( 788),P( 789)/  -79.07270,  -79.75015,  -82.16769/
      DATA P( 790),P( 791),P( 792)/  -82.74732,  -84.59505,  -82.60622/
      DATA P( 793),P( 794),P( 795)/  -81.71070,  -79.35495,  -77.74792/
      DATA P( 796),P( 797),P( 798)/  -74.77526,  -72.62601,  -68.55112/
      DATA P( 799),P( 800),P( 801)/  -65.83868,  -61.21409,  -58.36474/
      DATA P( 802),P( 803),P( 804)/  -54.30278,  -50.84036,  -46.69600/
      DATA P( 805),P( 806),P( 807)/  -43.59764,  -37.99600,  -31.69900/
      DATA P( 808),P( 809),P( 810)/  -40.69700,  -46.64900,  -54.39000/
      DATA P( 811),P( 812),P( 813)/  -57.97477,  -63.17451,  -65.47743/
      DATA P( 814),P( 815),P( 816)/  -70.30489,  -71.52654,  -76.00873/
      DATA P( 817),P( 818),P( 819)/  -76.79699,  -80.64429,  -81.10267/
      DATA P( 820),P( 821),P( 822)/  -84.52157,  -84.87836,  -87.91967/
      DATA P( 823),P( 824),P( 825)/  -86.20705,  -85.94186,  -83.63898/
      DATA P( 826),P( 827),P( 828)/  -82.87511,  -80.08761,  -78.84177/
      DATA P( 829),P( 830),P( 831)/  -75.11733,  -72.95416,  -68.79198/
      DATA P( 832),P( 833),P( 834)/  -66.62866,  -62.11663,  -60.21947/
      DATA P( 835),P( 836),P( 837)/  -55.40765,  -53.07765,  -47.55300/
      DATA P( 838),P( 839),P( 840)/  -44.40400,  -46.92900,  -52.62900/
      DATA P( 841),P( 842),P( 843)/  -58.35743,  -61.16500,  -66.01617/
      DATA P( 844),P( 845),P( 846)/  -68.19274,  -72.32810,  -74.15831/
      DATA P( 847),P( 848),P( 849)/  -77.84767,  -79.28156,  -83.01675/
      DATA P( 850),P( 851),P( 852)/  -84.29707,  -87.70210,  -86.48786/
      DATA P( 853),P( 854),P( 855)/  -86.34630,  -84.81548,  -84.22421/
      DATA P( 856),P( 857),P( 858)/  -82.34965,  -81.20419,  -78.34070/
      DATA P( 859),P( 860),P( 861)/  -76.26046,  -72.45203,  -70.20229/
      DATA P( 862),P( 863),P( 864)/  -67.29447,  -64.91265,  -61.89265/
      DATA P( 865),P( 866),P( 867)/  -58.74000,  -54.53900,  -51.14800/
      DATA P( 868),P( 869),P( 870)/  -46.37000,  -47.35700,  -55.37700/
      DATA P( 871),P( 872),P( 873)/  -58.85617,  -64.19273,  -66.46011/
      DATA P( 874),P( 875),P( 876)/  -71.49200,  -73.15467,  -77.80503/
      DATA P( 877),P( 878),P( 879)/  -79.34784,  -83.62377,  -84.86942/
      DATA P( 880),P( 881),P( 882)/  -88.76794,  -87.89114,  -88.45457/
      DATA P( 883),P( 884),P( 885)/  -87.11739,  -87.26630,  -85.65763/
      DATA P( 886),P( 887),P( 888)/  -85.44065,  -82.94886,  -81.27623/
      DATA P( 889),P( 890),P( 891)/  -77.76942,  -76.60447,  -73.45765/
      DATA P( 892),P( 893),P( 894)/  -71.74265,  -68.37439,  -66.34100/
      DATA P( 895),P( 896),P( 897)/  -62.36400,  -59.69900,  -55.08900/
      DATA P( 898),P( 899),P( 900)/  -51.90300,  -47.46000,  -52.97400/
      DATA P( 901),P( 902),P( 903)/  -58.96011,  -61.87900,  -67.15467/
      DATA P( 904),P( 905),P( 906)/  -69.82703,  -74.18284,  -76.42400/
      DATA P( 907),P( 908),P( 909)/  -80.57841,  -82.65694,  -86.63775/
      DATA P( 910),P( 911),P( 912)/  -86.44896,  -87.20875,  -86.36489/
      DATA P( 913),P( 914),P( 915)/  -86.78246,  -85.60423,  -85.60695/
      DATA P( 916),P( 917),P( 918)/  -83.52642,  -82.32742,  -79.93948/
      DATA P( 919),P( 920),P( 921)/  -78.94265,  -76.34764,  -75.31939/
      DATA P( 922),P( 923),P( 924)/  -72.22853,  -70.85500,  -66.89101/
      DATA P( 925),P( 926),P( 927)/  -64.91600,  -60.53800,  -58.09700/
      DATA P( 928),P( 929),P( 930)/  -53.39300,  -47.74800,  -55.80600/
      DATA P( 931),P( 932),P( 933)/  -59.06600,  -64.55703,  -67.69461/
      DATA P( 934),P( 935),P( 936)/  -72.70056,  -75.00336,  -80.16794/
      DATA P( 937),P( 938),P( 939)/  -82.20364,  -86.80547,  -86.80386/
      DATA P( 940),P( 941),P( 942)/  -88.41034,  -87.70808,  -88.79102/
      DATA P( 943),P( 944),P( 945)/  -87.54083,  -88.11201,  -85.96608/
      DATA P( 946),P( 947),P( 948)/  -86.18447,  -83.51165,  -83.55765/
      DATA P( 949),P( 950),P( 951)/  -80.84939,  -80.33353,  -77.34000/
      DATA P( 952),P( 953),P( 954)/  -76.25742,  -72.94089,  -71.18501/
      DATA P( 955),P( 956),P( 957)/  -67.24500,  -65.45600,  -61.00400/
      DATA P( 958),P( 959),P( 960)/  -58.83300,  -53.99900,  -47.56200/
      DATA P( 961),P( 962),P( 963)/  -53.20700,  -59.12200,  -62.56800/
      DATA P( 964),P( 965),P( 966)/  -67.49336,  -71.20728,  -75.98363/
      DATA P( 967),P( 968),P( 969)/  -78.93511,  -83.60300,  -84.15459/
      DATA P( 970),P( 971),P( 972)/  -86.01745,  -85.81779,  -87.22057/
      DATA P( 973),P( 974),P( 975)/  -86.42802,  -87.32331,  -86.01639/
      DATA P( 976),P( 977),P( 978)/  -86.33609,  -84.56760,  -84.59939/
      DATA P( 979),P( 980),P( 981)/  -82.48853,  -82.28999,  -79.77742/
      DATA P( 982),P( 983),P( 984)/  -79.10088,  -75.93541,  -74.86700/
      DATA P( 985),P( 986),P( 987)/  -71.36201,  -69.81500,  -65.91300/
      DATA P( 988),P( 989),P( 990)/  -63.96600,  -59.72700,  -57.49200/
      DATA P( 991),P( 992),P( 993)/  -47.33900,  -55.49800,  -59.51300/
      DATA P( 994),P( 995),P( 996)/  -65.40900,  -68.57897,  -74.40800/
      DATA P( 997),P( 998),P( 999)/  -77.26601,  -82.56802,  -83.45000/
      DATA P(1000),P(1001),P(1002)/  -86.07220,  -86.11238,  -88.22448/
      DATA P(1003),P(1004),P(1005)/  -87.61697,  -89.21880,  -87.94959/
      DATA P(1006),P(1007),P(1008)/  -89.09785,  -87.25893,  -88.09125/
      DATA P(1009),P(1010),P(1011)/  -85.92999,  -86.32442,  -83.92088/
      DATA P(1012),P(1013),P(1014)/  -83.65541,  -80.85221,  -80.13997/
      DATA P(1015),P(1016),P(1017)/  -76.79201,  -75.86700,  -72.15401/
      DATA P(1018),P(1019),P(1020)/  -70.79401,  -66.77901,  -65.05600/
      DATA P(1021),P(1022),P(1023)/  -60.74300,  -58.65600,  -47.15200/
      DATA P(1024),P(1025),P(1026)/  -53.21600,  -59.10300,  -63.36000/
      DATA P(1027),P(1028),P(1029)/  -69.17300,  -72.93800,  -78.34000/
      DATA P(1030),P(1031),P(1032)/  -79.62576,  -82.58938,  -83.16710/
      DATA P(1033),P(1034),P(1035)/  -85.57439,  -85.58880,  -87.40810/
      DATA P(1036),P(1037),P(1038)/  -86.77532,  -88.02228,  -86.95000/
      DATA P(1039),P(1040),P(1041)/  -87.84691,  -86.36382,  -86.86130/
      DATA P(1042),P(1043),P(1044)/  -85.01673,  -85.01221,  -82.94997/
      DATA P(1045),P(1046),P(1047)/  -82.28800,  -79.53700,  -78.78600/
      DATA P(1048),P(1049),P(1050)/  -75.59400,  -74.40338,  -71.06100/
      DATA P(1051),P(1052),P(1053)/  -69.53600,  -65.73600,  -63.93800/
      DATA P(1054),P(1055),P(1056)/  -59.82100,  -57.67800,  -47.05900/
      DATA P(1057),P(1058),P(1059)/  -55.49800,  -59.69900,  -66.35001/
      DATA P(1060),P(1061),P(1062)/  -70.15100,  -76.17576,  -77.79938/
      DATA P(1063),P(1064),P(1065)/  -81.30009,  -82.18780,  -85.22741/
      DATA P(1066),P(1067),P(1068)/  -85.42811,  -87.92584,  -87.47920/
      DATA P(1069),P(1070),P(1071)/  -89.39089,  -88.41363,  -89.90491/
      DATA P(1072),P(1073),P(1074)/  -88.37227,  -89.52173,  -87.60371/
      DATA P(1075),P(1076),P(1077)/  -88.34998,  -86.02910,  -86.33711/
      DATA P(1078),P(1079),P(1080)/  -83.69347,  -83.49416,  -80.40338/
      DATA P(1081),P(1082),P(1083)/  -79.96103,  -76.53201,  -75.46599/
      DATA P(1084),P(1085),P(1086)/  -72.02300,  -70.76601,  -66.90000/
      DATA P(1087),P(1088),P(1089)/  -65.39101,  -61.23600,  -53.30000/
      DATA P(1090),P(1091),P(1092)/  -60.10000,  -64.57101,  -70.79401/
      DATA P(1093),P(1094),P(1095)/  -72.88009,  -76.75780,  -78.18085/
      DATA P(1096),P(1097),P(1098)/  -81.22427,  -81.97147,  -84.79160/
      DATA P(1099),P(1100),P(1101)/  -85.11224,  -87.06838,  -86.93965/
      DATA P(1102),P(1103),P(1104)/  -88.40527,  -87.60355,  -88.71966/
      DATA P(1105),P(1106),P(1107)/  -87.45753,  -88.21742,  -86.62508/
      DATA P(1108),P(1109),P(1110)/  -87.03347,  -84.94488,  -84.98738/
      DATA P(1111),P(1112),P(1113)/  -82.56802,  -82.26565,  -79.56599/
      DATA P(1114),P(1115),P(1116)/  -78.55656,  -75.64792,  -74.65824/
      DATA P(1117),P(1118),P(1119)/  -71.42700,  -69.95500,  -66.57401/
      DATA P(1120),P(1121),P(1122)/  -64.70200,  -61.01300,  -58.79600/
      DATA P(1123),P(1124),P(1125)/  -56.10400,  -60.60300,  -67.46001/
      DATA P(1126),P(1127),P(1128)/  -69.85300,  -74.30505,  -75.74774/
      DATA P(1129),P(1130),P(1131)/  -79.38447,  -80.64971,  -83.97596/
      DATA P(1132),P(1133),P(1134)/  -84.33018,  -87.13380,  -86.98827/
      DATA P(1135),P(1136),P(1137)/  -89.25257,  -88.50536,  -90.34972/
      DATA P(1138),P(1139),P(1140)/  -89.25423,  -90.58105,  -89.04993/
      DATA P(1141),P(1142),P(1143)/  -90.02132,  -88.09087,  -88.71973/
      DATA P(1144),P(1145),P(1146)/  -86.42564,  -86.70891,  -83.90656/
      DATA P(1147),P(1148),P(1149)/  -83.97292,  -81.05824,  -80.57401/
      DATA P(1150),P(1151),P(1152)/  -77.31058,  -76.71011,  -73.35778/
      DATA P(1153),P(1154),P(1155)/  -72.32677,  -68.52552,  -67.29054/
      DATA P(1156),P(1157),P(1158)/  -63.09900,  -61.49700,  -53.80300/
      DATA P(1159),P(1160),P(1161)/  -60.91000,  -64.13425,  -68.40900/
      DATA P(1162),P(1163),P(1164)/  -70.13447,  -74.59972,  -76.06727/
      DATA P(1165),P(1166),P(1167)/  -79.48118,  -80.61043,  -83.56227/
      DATA P(1168),P(1169),P(1170)/  -84.09557,  -86.48541,  -86.47172/
      DATA P(1171),P(1172),P(1173)/  -88.38882,  -87.99512,  -89.36639/
      DATA P(1174),P(1175),P(1176)/  -88.56946,  -89.53675,  -88.24974/
      DATA P(1177),P(1178),P(1179)/  -88.94301,  -87.23009,  -87.70357/
      DATA P(1180),P(1181),P(1182)/  -85.73330,  -85.83824,  -83.57633/
      DATA P(1183),P(1184),P(1185)/  -83.42558,  -80.87611,  -80.47977/
      DATA P(1186),P(1187),P(1188)/  -77.81277,  -76.99352,  -74.36054/
      DATA P(1189),P(1190),P(1191)/  -72.97514,  -69.99717,  -68.21571/
      DATA P(1192),P(1193),P(1194)/  -62.48554,  -57.43600,  -51.54900/
      DATA P(1195),P(1196),P(1197)/  -56.86400,  -59.56000,  -64.74800/
      DATA P(1198),P(1199),P(1200)/  -66.94601,  -71.55227,  -73.22436/
      DATA P(1201),P(1202),P(1203)/  -77.42533,  -78.55595,  -82.00375/
      DATA P(1204),P(1205),P(1206)/  -82.63573,  -85.83466,  -85.94391/
      DATA P(1207),P(1208),P(1209)/  -88.65883,  -88.33043,  -90.55814/
      DATA P(1210),P(1211),P(1212)/  -90.03264,  -91.52473,  -90.39797/
      DATA P(1213),P(1214),P(1215)/  -91.65311,  -90.06719,  -91.10330/
      DATA P(1216),P(1217),P(1218)/  -89.20277,  -89.94492,  -87.81947/
      DATA P(1219),P(1220),P(1221)/  -88.23611,  -85.89777,  -86.01978/
      DATA P(1222),P(1223),P(1224)/  -83.50752,  -83.33614,  -80.63013/
      DATA P(1225),P(1226),P(1227)/  -80.24617,  -77.38931,  -76.62054/
      DATA P(1228),P(1229),P(1230)/  -70.96671,  -66.63574,  -60.79900/
      DATA P(1231),P(1232),P(1233)/  -56.50400,  -50.49600,  -55.77800/
      DATA P(1234),P(1235),P(1236)/  -59.34800,  -63.78070,  -66.35700/
      DATA P(1237),P(1238),P(1239)/  -70.65401,  -72.50700,  -76.25573/
      DATA P(1240),P(1241),P(1242)/  -77.53500,  -80.84400,  -81.60386/
      DATA P(1243),P(1244),P(1245)/  -84.41389,  -84.67664,  -87.00264/
      DATA P(1246),P(1247),P(1248)/  -86.81781,  -88.64134,  -87.99648/
      DATA P(1249),P(1250),P(1251)/  -89.47328,  -88.42270,  -89.59290/
      DATA P(1252),P(1253),P(1254)/  -88.32852,  -89.22250,  -87.61862/
      DATA P(1255),P(1256),P(1257)/  -88.26109,  -86.39777,  -86.70852/
      DATA P(1258),P(1259),P(1260)/  -84.61008,  -84.62614,  -82.39394/
      DATA P(1261),P(1262),P(1263)/  -82.02131,  -79.72353,  -78.95672/
      DATA P(1264),P(1265),P(1266)/  -74.00574,  -69.70559,  -64.59000/
      DATA P(1267),P(1268),P(1269)/  -60.25800,  -54.99500,  -50.57100/
      DATA P(1270),P(1271),P(1272)/  -58.03000,  -60.51300,  -65.68258/
      DATA P(1273),P(1274),P(1275)/  -67.57384,  -72.27725,  -73.47569/
      DATA P(1276),P(1277),P(1278)/  -77.25660,  -78.31400,  -81.91901/
      DATA P(1279),P(1280),P(1281)/  -82.36398,  -85.30597,  -85.10670/
      DATA P(1282),P(1283),P(1284)/  -87.72326,  -87.18027,  -89.40488/
      DATA P(1285),P(1286),P(1287)/  -88.55730,  -90.31107,  -89.16916/
      DATA P(1288),P(1289),P(1290)/  -90.52307,  -89.02779,  -90.07030/
      DATA P(1291),P(1292),P(1293)/  -88.28952,  -88.99364,  -87.00564/
      DATA P(1294),P(1295),P(1296)/  -87.35294,  -85.21131,  -85.20953/
      DATA P(1297),P(1298),P(1299)/  -82.95972,  -82.39944,  -77.82559/
      DATA P(1300),P(1301),P(1302)/  -74.42342,  -69.55952,  -65.93100/
      DATA P(1303),P(1304),P(1305)/  -60.79900,  -57.10100,  -51.80000/
      DATA P(1306),P(1307),P(1308)/  -47.97200,  -52.82400,  -57.57409/
      DATA P(1309),P(1310),P(1311)/  -60.34800,  -64.94701,  -67.09600/
      DATA P(1312),P(1313),P(1314)/  -71.12492,  -72.79601,  -76.45901/
      DATA P(1315),P(1316),P(1317)/  -77.56082,  -80.43665,  -80.69045/
      DATA P(1318),P(1319),P(1320)/  -83.66600,  -83.78989,  -86.28795/
      DATA P(1321),P(1322),P(1323)/  -86.07706,  -87.93494,  -87.36349/
      DATA P(1324),P(1325),P(1326)/  -88.84202,  -87.91496,  -88.98708/
      DATA P(1327),P(1328),P(1329)/  -87.74183,  -88.50358,  -86.93258/
      DATA P(1330),P(1331),P(1332)/  -87.44476,  -85.70254,  -85.87772/
      DATA P(1333),P(1334),P(1335)/  -83.94944,  -83.78759,  -79.49825/
      DATA P(1336),P(1337),P(1338)/  -76.50112,  -72.29914,  -68.84354/
      DATA P(1339),P(1340),P(1341)/  -64.07700,  -60.70500,  -55.72200/
      DATA P(1342),P(1343),P(1344)/  -52.09800,  -46.93800,  -51.72100/
      DATA P(1345),P(1346),P(1347)/  -54.36900,  -59.92733,  -62.05348/
      DATA P(1348),P(1349),P(1350)/  -66.93201,  -68.43000,  -72.90100/
      DATA P(1351),P(1352),P(1353)/  -73.99381,  -77.71368,  -78.66066/
      DATA P(1354),P(1355),P(1356)/  -81.82989,  -82.54294,  -85.18661/
      DATA P(1357),P(1358),P(1359)/  -85.25894,  -87.65751,  -87.18948/
      DATA P(1360),P(1361),P(1362)/  -89.17297,  -88.32465,  -89.86081/
      DATA P(1363),P(1364),P(1365)/  -88.69736,  -89.88180,  -88.41561/
      DATA P(1366),P(1367),P(1368)/  -89.27953,  -87.64830,  -88.12444/
      DATA P(1369),P(1370),P(1371)/  -86.43565,  -86.42445,  -82.37859/
      DATA P(1372),P(1373),P(1374)/  -80.11914,  -75.64954,  -72.99590/
      DATA P(1375),P(1376),P(1377)/  -68.32854,  -65.48125,  -60.65000/
      DATA P(1378),P(1379),P(1380)/  -57.53800,  -52.47100,  -49.09000/
      DATA P(1381),P(1382),P(1383)/  -43.77100,  -46.26600,  -51.66483/
      DATA P(1384),P(1385),P(1386)/  -54.56600,  -59.67200,  -62.49006/
      DATA P(1387),P(1388),P(1389)/  -66.47189,  -68.41368,  -72.31105/
      DATA P(1390),P(1391),P(1392)/  -73.88776,  -77.14294,  -78.13190/
      DATA P(1393),P(1394),P(1395)/  -81.04913,  -81.74257,  -84.09073/
      DATA P(1396),P(1397),P(1398)/  -84.34869,  -86.23994,  -85.93225/
      DATA P(1399),P(1400),P(1401)/  -87.50140,  -86.90263,  -88.06322/
      DATA P(1402),P(1403),P(1404)/  -87.16007,  -88.07566,  -86.89588/
      DATA P(1405),P(1406),P(1407)/  -87.58660,  -86.34414,  -86.55116/
      DATA P(1408),P(1409),P(1410)/  -82.89314,  -80.70657,  -77.05590/
      DATA P(1411),P(1412),P(1413)/  -74.47855,  -70.52125,  -67.69139/
      DATA P(1414),P(1415),P(1416)/  -63.31609,  -60.18547,  -55.73870/
      DATA P(1417),P(1418),P(1419)/  -52.28994,  -47.59977,  -44.04100/
      DATA P(1420),P(1421),P(1422)/  -39.15100,  -35.39700,  -45.69800/
      DATA P(1423),P(1424),P(1425)/  -48.70800,  -54.32500,  -56.95200/
      DATA P(1426),P(1427),P(1428)/  -62.00000,  -64.22471,  -68.88776/
      DATA P(1429),P(1430),P(1431)/  -70.34092,  -74.27700,  -75.59100/
      DATA P(1432),P(1433),P(1434)/  -79.09460,  -79.53073,  -82.67554/
      DATA P(1435),P(1436),P(1437)/  -82.78994,  -85.40973,  -85.06985/
      DATA P(1438),P(1439),P(1440)/  -87.27122,  -86.69340,  -88.43961/
      DATA P(1441),P(1442),P(1443)/  -87.55822,  -88.95455,  -87.85594/
      DATA P(1444),P(1445),P(1446)/  -88.89236,  -87.72678,  -88.26717/
      DATA P(1447),P(1448),P(1449)/  -84.91929,  -83.27603,  -79.72988/
      DATA P(1450),P(1451),P(1452)/  -77.82802,  -73.94461,  -71.78049/
      DATA P(1453),P(1454),P(1455)/  -68.07003,  -65.10523,  -61.48557/
      DATA P(1456),P(1457),P(1458)/  -58.04849,  -53.59800,  -50.65500/
      DATA P(1459),P(1460),P(1461)/  -45.92300,  -42.70000,  -37.62300/
      DATA P(1462),P(1463),P(1464)/  -46.56500,  -49.77000,  -54.96700/
      DATA P(1465),P(1466),P(1467)/  -57.68700,  -62.40100,  -64.54300/
      DATA P(1468),P(1469),P(1470)/  -68.70700,  -70.30000,  -73.89500/
      DATA P(1471),P(1472),P(1473)/  -75.10600,  -78.09600,  -78.75973/
      DATA P(1474),P(1475),P(1476)/  -81.34985,  -81.67300,  -83.73339/
      DATA P(1477),P(1478),P(1479)/  -83.73161,  -85.32822,  -85.24137/
      DATA P(1480),P(1481),P(1482)/  -86.65594,  -86.02236,  -87.12667/
      DATA P(1483),P(1484),P(1485)/  -86.52943,  -87.23612,  -84.32577/
      DATA P(1486),P(1487),P(1488)/  -82.94299,  -80.03909,  -78.19086/
      DATA P(1489),P(1490),P(1491)/  -74.89987,  -72.99338,  -69.20986/
      DATA P(1492),P(1493),P(1494)/  -67.23557,  -63.16349,  -61.13400/
      DATA P(1495),P(1496),P(1497)/  -57.22200,  -54.43700,  -50.19800/
      DATA P(1498),P(1499),P(1500)/  -47.08700,  -42.47600,  -39.00200/
      DATA P(1501),P(1502),P(1503)/  -44.00400,  -49.70500,  -52.47100/
      DATA P(1504),P(1505),P(1506)/  -57.74300,  -60.07200,  -64.72000/
      DATA P(1507),P(1508),P(1509)/  -66.56500,  -70.70000,  -71.95800/
      DATA P(1510),P(1511),P(1512)/  -75.57201,  -76.30000,  -79.46201/
      DATA P(1513),P(1514),P(1515)/  -79.71339,  -82.44701,  -82.39101/
      DATA P(1516),P(1517),P(1518)/  -84.74137,  -84.63036,  -86.49519/
      DATA P(1519),P(1520),P(1521)/  -85.90457,  -87.57386,  -86.95812/
      DATA P(1522),P(1523),P(1524)/  -88.08762,  -85.44491,  -84.54263/
      DATA P(1525),P(1526),P(1527)/  -81.61642,  -80.44131,  -77.10173/
      DATA P(1528),P(1529),P(1530)/  -75.74003,  -72.18056,  -70.42584/
      DATA P(1531),P(1532),P(1533)/  -66.79816,  -64.99368,  -61.44100/
      DATA P(1534),P(1535),P(1536)/  -59.26200,  -55.34900,  -52.79700/
      DATA P(1537),P(1538),P(1539)/  -48.40000,  -45.40100,  -40.66900/
      DATA P(1540),P(1541),P(1542)/  -41.57900,  -45.03800,  -50.33800/
      DATA P(1543),P(1544),P(1545)/  -53.13200,  -57.91100,  -60.25800/
      DATA P(1546),P(1547),P(1548)/  -64.43100,  -66.32201,  -69.99200/
      DATA P(1549),P(1550),P(1551)/  -71.37100,  -74.46339,  -75.33900/
      DATA P(1552),P(1553),P(1554)/  -78.05901,  -78.55100,  -80.91036/
      DATA P(1555),P(1556),P(1557)/  -81.36885,  -83.20258,  -83.13686/
      DATA P(1558),P(1559),P(1560)/  -84.82912,  -84.69962,  -86.02558/
      DATA P(1561),P(1562),P(1563)/  -83.79732,  -83.07787,  -80.75996/
      DATA P(1564),P(1565),P(1566)/  -79.63631,  -76.76626,  -75.47057/
      DATA P(1567),P(1568),P(1569)/  -72.48584,  -70.98817,  -68.00368/
      DATA P(1570),P(1571),P(1572)/  -66.85530,  -63.71400,  -61.80500/
      DATA P(1573),P(1574),P(1575)/  -58.32100,  -55.89900,  -52.05200/
      DATA P(1576),P(1577),P(1578)/  -49.21100,  -44.91700,  -41.70300/
      DATA P(1579),P(1580),P(1581)/  -53.03000,  -55.42400,  -60.18400/
      DATA P(1582),P(1583),P(1584)/  -62.17300,  -66.34100,  -67.90340/
      DATA P(1585),P(1586),P(1587)/  -71.61301,  -72.46101,  -75.78101/
      DATA P(1588),P(1589),P(1590)/  -76.15900,  -79.15784,  -79.51257/
      DATA P(1591),P(1592),P(1593)/  -82.03700,  -82.04212,  -84.47735/
      DATA P(1594),P(1595),P(1596)/  -84.20258,  -85.95952,  -84.01179/
      DATA P(1597),P(1598),P(1599)/  -83.75748,  -81.44158,  -80.93551/
      DATA P(1600),P(1601),P(1602)/  -78.15625,  -77.41784,  -74.38520/
      DATA P(1603),P(1604),P(1605)/  -73.69368,  -70.95679,  -70.15786/
      DATA P(1606),P(1607),P(1608)/  -67.35210,  -65.68587,  -62.75516/
      DATA P(1609),P(1610),P(1611)/  -60.36100,  -56.57000,  -54.14800/
      DATA P(1612),P(1613),P(1614)/  -49.93700,  -47.14300,  -42.54100/
      DATA P(1615),P(1616),P(1617)/  -48.19500,  -52.94600,  -55.47000/
      DATA P(1618),P(1619),P(1620)/  -59.80200,  -61.71100,  -65.46500/
      DATA P(1621),P(1622),P(1623)/  -66.61100,  -70.21900,  -71.30785/
      DATA P(1624),P(1625),P(1626)/  -73.85600,  -75.03700,  -77.53772/
      DATA P(1627),P(1628),P(1629)/  -78.43025,  -80.47489,  -81.08585/
      DATA P(1630),P(1631),P(1632)/  -82.97043,  -81.42583,  -81.27854/
      DATA P(1633),P(1634),P(1635)/  -79.46373,  -79.05225,  -76.87825/
      DATA P(1636),P(1637),P(1638)/  -76.07585,  -73.60714,  -73.39921/
      DATA P(1639),P(1640),P(1641)/  -71.26808,  -70.68810,  -68.42101/
      DATA P(1642),P(1643),P(1644)/  -66.97716,  -64.21686,  -62.22400/
      DATA P(1645),P(1646),P(1647)/  -58.97300,  -56.70000,  -53.10400/
      DATA P(1648),P(1649),P(1650)/  -50.43100,  -46.30500,  -43.29600/
      DATA P(1651),P(1652),P(1653)/  -47.85100,  -50.40300,  -55.12600/
      DATA P(1654),P(1655),P(1656)/  -57.07300,  -61.46000,  -63.01600/
      DATA P(1657),P(1658),P(1659)/  -66.78800,  -67.95555,  -71.22200/
      DATA P(1660),P(1661),P(1662)/  -72.37521,  -75.45939,  -75.94608/
      DATA P(1663),P(1664),P(1665)/  -78.99695,  -79.52763,  -81.97637/
      DATA P(1666),P(1667),P(1668)/  -80.66215,  -81.00573,  -79.27640/
      DATA P(1669),P(1670),P(1671)/  -79.34660,  -77.14677,  -77.06113/
      DATA P(1672),P(1673),P(1674)/  -74.58625,  -74.77265,  -72.56905/
      DATA P(1675),P(1676),P(1677)/  -72.46529,  -70.20116,  -69.37186/
      DATA P(1678),P(1679),P(1680)/  -66.73734,  -65.21581,  -62.22400/
      DATA P(1681),P(1682),P(1683)/  -60.41700,  -56.97900,  -54.75300/
      DATA P(1684),P(1685),P(1686)/  -50.89700,  -48.17700,  -43.79900/
      DATA P(1687),P(1688),P(1689)/  -42.70000,  -47.59900,  -50.00300/
      DATA P(1690),P(1691),P(1692)/  -54.28700,  -56.35500,  -60.35100/
      DATA P(1693),P(1694),P(1695)/  -61.99100,  -65.35500,  -66.98939/
      DATA P(1696),P(1697),P(1698)/  -69.96835,  -71.35240,  -74.25251/
      DATA P(1699),P(1700),P(1701)/  -75.66142,  -78.00211,  -77.12797/
      DATA P(1702),P(1703),P(1704)/  -77.55506,  -76.23926,  -76.45150/
      DATA P(1705),P(1706),P(1707)/  -74.80055,  -74.66294,  -72.89834/
      DATA P(1708),P(1709),P(1710)/  -73.37730,  -71.74796,  -71.82803/
      DATA P(1711),P(1712),P(1713)/  -70.09412,  -69.47134,  -67.21481/
      DATA P(1714),P(1715),P(1716)/  -66.05736,  -63.37200,  -61.77700/
      DATA P(1717),P(1718),P(1719)/  -58.64700,  -56.62600,  -53.10400/
      DATA P(1720),P(1721),P(1722)/  -50.56100,  -46.60300,  -43.73400/
      DATA P(1723),P(1724),P(1725)/  -49.30400,  -51.55800,  -55.91800/
      DATA P(1726),P(1727),P(1728)/  -57.67800,  -61.52900,  -63.14600/
      DATA P(1729),P(1730),P(1731)/  -66.85201,  -68.24251,  -71.92101/
      DATA P(1732),P(1733),P(1734)/  -72.94762,  -76.09808,  -75.36769/
      DATA P(1735),P(1736),P(1737)/  -76.28024,  -75.13763,  -75.77195/
      DATA P(1738),P(1739),P(1740)/  -74.19882,  -74.71709,  -72.89286/
      DATA P(1741),P(1742),P(1743)/  -73.71631,  -72.08012,  -72.54516/
      DATA P(1744),P(1745),P(1746)/  -70.83389,  -70.69990,  -68.57185/
      DATA P(1747),P(1748),P(1749)/  -67.95191,  -65.51598,  -64.29058/
      DATA P(1750),P(1751),P(1752)/  -61.48800,  -59.74600,  -56.46700/
      DATA P(1753),P(1754),P(1755)/  -54.39900,  -50.70100,  -48.10200/
      DATA P(1756),P(1757),P(1758)/  -43.90100,  -43.90100,  -48.41000/
      DATA P(1759),P(1760),P(1761)/  -50.72900,  -54.80900,  -56.95200/
      DATA P(1762),P(1763),P(1764)/  -60.78000,  -62.84800,  -66.24800/
      DATA P(1765),P(1766),P(1767)/  -67.83080,  -70.75891,  -70.51537/
      DATA P(1768),P(1769),P(1770)/  -71.49995,  -71.11569,  -71.63359/
      DATA P(1771),P(1772),P(1773)/  -70.72710,  -71.32369,  -70.15430/
      DATA P(1774),P(1775),P(1776)/  -71.25890,  -70.10074,  -70.77383/
      DATA P(1777),P(1778),P(1779)/  -69.47989,  -69.54242,  -67.84627/
      DATA P(1780),P(1781),P(1782)/  -67.47156,  -65.68448,  -64.60474/
      DATA P(1783),P(1784),P(1785)/  -62.08663,  -60.65900,  -57.70600/
      DATA P(1786),P(1787),P(1788)/  -55.84300,  -52.49900,  -50.09600/
      DATA P(1789),P(1790),P(1791)/  -46.34200,  -43.50100,  -43.04400/
      DATA P(1792),P(1793),P(1794)/  -45.46600,  -50.05200,  -52.32200/
      DATA P(1795),P(1796),P(1797)/  -56.75600,  -58.72800,  -62.67080/
      DATA P(1798),P(1799),P(1800)/  -64.38626,  -67.83337,  -67.68796/
      DATA P(1801),P(1802),P(1803)/  -69.32203,  -68.76322,  -70.12856/
      DATA P(1804),P(1805),P(1806)/  -69.15330,  -70.40040,  -69.16440/
      DATA P(1807),P(1808),P(1809)/  -70.53432,  -69.43239,  -70.41662/
      DATA P(1810),P(1811),P(1812)/  -69.17677,  -69.68159,  -68.06464/
      DATA P(1813),P(1814),P(1815)/  -68.19026,  -66.38987,  -65.97663/
      DATA P(1816),P(1817),P(1818)/  -63.62119,  -62.59337,  -59.94254/
      DATA P(1819),P(1820),P(1821)/  -58.47000,  -55.60679,  -53.40300/
      DATA P(1822),P(1823),P(1824)/  -49.85400,  -47.40400,  -43.37000/
      DATA P(1825),P(1826),P(1827)/  -37.39000,  -42.20600,  -45.04700/
      DATA P(1828),P(1829),P(1830)/  -49.48100,  -52.07100,  -56.03900/
      DATA P(1831),P(1832),P(1833)/  -58.43300,  -61.67426,  -62.08200/
      DATA P(1834),P(1835),P(1836)/  -63.63892,  -63.58322,  -65.02339/
      DATA P(1837),P(1838),P(1839)/  -64.64915,  -66.06240,  -65.47401/
      DATA P(1840),P(1841),P(1842)/  -66.89239,  -66.18662,  -67.33906/
      DATA P(1843),P(1844),P(1845)/  -66.39159,  -67.20573,  -66.04998/
      DATA P(1846),P(1847),P(1848)/  -66.38731,  -64.98979,  -64.90727/
      DATA P(1849),P(1850),P(1851)/  -63.07959,  -62.29254,  -60.08469/
      DATA P(1852),P(1853),P(1854)/  -58.80679,  -56.24830,  -54.52851/
      DATA P(1855),P(1856),P(1857)/  -51.40000,  -49.09900,  -45.50300/
      DATA P(1858),P(1859),P(1860)/  -42.80200,  -36.71000,  -39.62600/
      DATA P(1861),P(1862),P(1863)/  -44.60000,  -47.21700,  -51.75400/
      DATA P(1864),P(1865),P(1866)/  -53.86400,  -57.97400,  -58.25600/
      DATA P(1867),P(1868),P(1869)/  -60.47403,  -60.46036,  -62.61754/
      DATA P(1870),P(1871),P(1872)/  -62.21981,  -64.25910,  -63.39233/
      DATA P(1873),P(1874),P(1875)/  -65.28700,  -64.57049,  -66.06255/
      DATA P(1876),P(1877),P(1878)/  -65.20332,  -66.34572,  -65.17731/
      DATA P(1879),P(1880),P(1881)/  -65.95257,  -64.53129,  -64.93447/
      DATA P(1882),P(1883),P(1884)/  -63.29925,  -62.99900,  -60.93081/
      DATA P(1885),P(1886),P(1887)/  -60.11831,  -57.72851,  -56.49310/
      DATA P(1888),P(1889),P(1890)/  -53.65400,  -51.84700,  -48.50300/
      DATA P(1891),P(1892),P(1893)/  -46.30500,  -42.50400,  -31.21000/
      DATA P(1894),P(1895),P(1896)/  -36.25300,  -39.54200,  -44.10600/
      DATA P(1897),P(1898),P(1899)/  -46.88200,  -50.82800,  -51.88400/
      DATA P(1900),P(1901),P(1902)/  -54.00091,  -54.56400,  -56.64269/
      DATA P(1903),P(1904),P(1905)/  -56.81470,  -58.91133,  -58.68700/
      DATA P(1906),P(1907),P(1908)/  -60.72505,  -60.46255,  -62.03932/
      DATA P(1909),P(1910),P(1911)/  -61.50640,  -62.73831,  -61.99001/
      DATA P(1912),P(1913),P(1914)/  -62.93875,  -61.89485,  -62.55089/
      DATA P(1915),P(1916),P(1917)/  -61.31989,  -61.28194,  -59.80389/
      DATA P(1918),P(1919),P(1920)/  -59.21896,  -57.38364,  -56.26192/
      DATA P(1921),P(1922),P(1923)/  -53.87331,  -52.31932,  -49.37718/
      DATA P(1924),P(1925),P(1926)/  -47.46900,  -44.11600,  -41.60100/
      DATA P(1927),P(1928),P(1929)/  -30.96300,  -34.01800,  -39.13200/
      DATA P(1930),P(1931),P(1932)/  -41.68500,  -46.41900,  -47.31100/
      DATA P(1933),P(1934),P(1935)/  -50.07500,  -50.49400,  -53.23757/
      DATA P(1936),P(1937),P(1938)/  -53.41312,  -56.02200,  -55.74643/
      DATA P(1939),P(1940),P(1941)/  -58.16300,  -57.88900,  -59.84800/
      DATA P(1942),P(1943),P(1944)/  -59.36831,  -60.99400,  -60.17675/
      DATA P(1945),P(1946),P(1947)/  -61.59073,  -60.59660,  -61.57690/
      DATA P(1948),P(1949),P(1950)/  -60.37280,  -60.77192,  -59.31539/
      DATA P(1951),P(1952),P(1953)/  -59.26379,  -57.56003,  -56.95331/
      DATA P(1954),P(1955),P(1956)/  -54.70432,  -53.49718,  -50.99265/
      DATA P(1957),P(1958),P(1959)/  -49.70135,  -46.41600,  -44.40400/
      DATA P(1960),P(1961),P(1962)/  -40.84600,  -25.46000,  -30.60200/
      DATA P(1963),P(1964),P(1965)/  -33.89700,  -38.48000,  -39.96100/
      DATA P(1966),P(1967),P(1968)/  -42.63200,  -43.86600,  -46.48011/
      DATA P(1969),P(1970),P(1971)/  -47.34900,  -49.72770,  -50.28300/
      DATA P(1972),P(1973),P(1974)/  -52.58900,  -52.88800,  -54.76831/
      DATA P(1975),P(1976),P(1977)/  -54.75800,  -56.25675,  -56.11073/
      DATA P(1978),P(1979),P(1980)/  -57.46660,  -57.10190,  -58.07980/
      DATA P(1981),P(1982),P(1983)/  -57.31278,  -57.83654,  -56.74452/
      DATA P(1984),P(1985),P(1986)/  -56.88922,  -55.57896,  -55.17434/
      DATA P(1987),P(1988),P(1989)/  -53.39100,  -52.39188,  -50.34597/
      DATA P(1990),P(1991),P(1992)/  -49.06717,  -46.68650,  -44.74000/
      DATA P(1993),P(1994),P(1995)/  -41.72200,  -39.52300,  -36.17000/
      DATA P(1996),P(1997),P(1998)/  -33.30100,  -34.68900,  -37.96100/
      DATA P(1999),P(2000),P(2001)/  -39.00400,  -42.24600,  -42.84600/
      DATA P(2002),P(2003),P(2004)/  -45.90999,  -46.26651,  -49.18010/
      DATA P(2005),P(2006),P(2007)/  -49.31600,  -51.77000,  -51.66100/
      DATA P(2008),P(2009),P(2010)/  -53.79400,  -53.46800,  -55.30300/
      DATA P(2011),P(2012),P(2013)/  -54.81044,  -56.21600,  -55.43300/
      DATA P(2014),P(2015),P(2016)/  -56.39452,  -55.28400,  -55.85223/
      DATA P(2017),P(2018),P(2019)/  -54.48961,  -54.58384,  -52.89021/
      DATA P(2020),P(2021),P(2022)/  -52.44522,  -50.47293,  -49.78951/
      DATA P(2023),P(2024),P(2025)/  -47.41385,  -46.05968,  -43.28558/
      DATA P(2026),P(2027),P(2028)/  -41.50003,  -38.39600,  -36.40300/
      DATA P(2029),P(2030),P(2031)/  -26.37100,  -29.67300,  -31.32800/
      DATA P(2032),P(2033),P(2034)/  -34.54500,  -35.99500,  -38.77530/
      DATA P(2035),P(2036),P(2037)/  -39.91700,  -42.55376,  -43.24900/
      DATA P(2038),P(2039),P(2040)/  -45.81300,  -46.13700,  -48.46300/
      DATA P(2041),P(2042),P(2043)/  -48.63600,  -50.37500,  -50.21700/
      DATA P(2044),P(2045),P(2046)/  -51.73500,  -51.47452,  -52.59400/
      DATA P(2047),P(2048),P(2049)/  -52.00723,  -52.49000,  -51.47384/
      DATA P(2050),P(2051),P(2052)/  -51.72421,  -50.53322,  -50.36204/
      DATA P(2053),P(2054),P(2055)/  -48.93542,  -48.44109,  -46.43272/
      DATA P(2056),P(2057),P(2058)/  -45.29558,  -42.84003,  -41.39644/
      DATA P(2059),P(2060),P(2061)/  -38.61033,  -36.87800,  -33.80400/
      DATA P(2062),P(2063),P(2064)/  -24.27600,  -25.82100,  -29.46400/
      DATA P(2065),P(2066),P(2067)/  -30.65600,  -34.14700,  -34.90100/
      DATA P(2068),P(2069),P(2070)/  -38.20633,  -38.80980,  -41.89869/
      DATA P(2071),P(2072),P(2073)/  -42.22300,  -44.83900,  -44.93600/
      DATA P(2074),P(2075),P(2076)/  -47.23600,  -47.07800,  -48.97500/
      DATA P(2077),P(2078),P(2079)/  -48.59400,  -50.15200,  -49.58300/
      DATA P(2080),P(2081),P(2082)/  -50.68300,  -49.72300,  -50.44192/
      DATA P(2083),P(2084),P(2085)/  -49.30236,  -49.64328,  -48.25320/
      DATA P(2086),P(2087),P(2088)/  -48.24625,  -46.36561,  -45.70603/
      DATA P(2089),P(2090),P(2091)/  -43.38844,  -42.51133,  -39.90672/
      DATA P(2092),P(2093),P(2094)/  -38.66915,  -35.47854,  -34.29803/
      DATA P(2095),P(2096),P(2097)/  -17.24700,  -20.80900,  -22.62900/
      DATA P(2098),P(2099),P(2100)/  -26.11200,  -27.64700,  -30.69247/
      DATA P(2101),P(2102),P(2103)/  -31.85500,  -34.87200,  -35.76100/
      DATA P(2104),P(2105),P(2106)/  -38.35000,  -38.97100,  -41.40800/
      DATA P(2107),P(2108),P(2109)/  -41.65100,  -43.72200,  -43.67600/
      DATA P(2110),P(2111),P(2112)/  -45.27700,  -45.11200,  -46.32300/
      DATA P(2113),P(2114),P(2115)/  -45.78192,  -46.59236,  -45.84098/
      DATA P(2116),P(2117),P(2118)/  -46.51453,  -45.44624,  -45.80961/
      DATA P(2119),P(2120),P(2121)/  -44.22334,  -43.82143,  -41.92977/
      DATA P(2122),P(2123),P(2124)/  -41.21787,  -39.01815,  -37.97854/
      DATA P(2125),P(2126),P(2127)/  -35.56803,  -34.35015,  -31.70800/
      DATA P(2128),P(2129),P(2130)/  -15.07200,  -16.72200,  -20.56100/
      DATA P(2131),P(2132),P(2133)/  -21.91400,  -25.59200,  -26.49700/
      DATA P(2134),P(2135),P(2136)/  -29.96345,  -30.66831,  -33.93459/
      DATA P(2137),P(2138),P(2139)/  -34.42800,  -37.18700,  -37.45400/
      DATA P(2140),P(2141),P(2142)/  -39.93900,  -39.98000,  -41.96400/
      DATA P(2143),P(2144),P(2145)/  -41.87500,  -43.45584,  -42.89400/
      DATA P(2146),P(2147),P(2148)/  -44.38500,  -43.52452,  -44.53825/
      DATA P(2149),P(2150),P(2151)/  -43.67800,  -44.25452,  -42.80864/
      DATA P(2152),P(2153),P(2154)/  -42.99929,  -41.22054,  -41.13850/
      DATA P(2155),P(2156),P(2157)/  -38.98780,  -38.70803,  -36.39538/
      DATA P(2158),P(2159),P(2160)/  -35.88203,  -33.39584,  -32.43526/
      DATA P(2161),P(2162),P(2163)/  -29.69240,  -28.29641,  -11.56900/
      DATA P(2164),P(2165),P(2166)/  -13.50100,  -17.19300,  -18.66200/
      DATA P(2167),P(2168),P(2169)/  -21.99176,  -23.25700,  -26.28800/
      DATA P(2170),P(2171),P(2172)/  -27.34600,  -30.08000,  -30.92200/
      DATA P(2173),P(2174),P(2175)/  -33.27400,  -33.98900,  -36.17000/
      DATA P(2176),P(2177),P(2178)/  -36.25100,  -38.05200,  -37.95500/
      DATA P(2179),P(2180),P(2181)/  -39.45607,  -39.00380,  -40.22800/
      DATA P(2182),P(2183),P(2184)/  -39.69252,  -40.43600,  -39.16829/
      DATA P(2185),P(2186),P(2187)/  -39.71813,  -38.32915,  -38.45535/
      DATA P(2188),P(2189),P(2190)/  -36.70803,  -36.70906,  -34.83583/
      DATA P(2191),P(2192),P(2193)/  -34.53635,  -32.53186,  -31.69240/
      DATA P(2194),P(2195),P(2196)/  -29.45393,  -28.28342,  -25.82100/
      DATA P(2197),P(2198),P(2199)/  -24.41710,  -11.14600,  -12.64900/
      DATA P(2200),P(2201),P(2202)/  -16.46300,  -17.46500,  -21.07399/
      DATA P(2203),P(2204),P(2205)/  -21.89044,  -25.32613,  -25.82500/
      DATA P(2206),P(2207),P(2208)/  -28.87600,  -29.38600,  -31.94100/
      DATA P(2209),P(2210),P(2211)/  -32.16000,  -34.26500,  -34.30000/
      DATA P(2212),P(2213),P(2214)/  -36.07896,  -35.65000,  -37.35800/
      DATA P(2215),P(2216),P(2217)/  -36.55761,  -37.78852,  -36.74000/
      DATA P(2218),P(2219),P(2220)/  -37.82266,  -36.48484,  -37.32489/
      DATA P(2221),P(2222),P(2223)/  -35.69051,  -36.29551,  -34.47971/
      DATA P(2224),P(2225),P(2226)/  -34.77864,  -32.81239,  -32.66294/
      DATA P(2227),P(2228),P(2229)/  -30.43805,  -29.92330,  -27.40808/
      DATA P(2230),P(2231),P(2232)/  -26.61847,  -23.75639,  -22.59800/
      DATA P(2233),P(2234),P(2235)/   -7.66200,   -9.21300,  -12.67005/
      DATA P(2236),P(2237),P(2238)/  -14.05000,  -17.18500,  -18.37900/
      DATA P(2239),P(2240),P(2241)/  -21.22400,  -22.37900,  -24.76700/
      DATA P(2242),P(2243),P(2244)/  -25.71300,  -27.99300,  -28.29900/
      DATA P(2245),P(2246),P(2247)/  -30.16100,  -30.29800,  -31.85061/
      DATA P(2248),P(2249),P(2250)/  -31.67296,  -33.01000,  -32.52300/
      DATA P(2251),P(2252),P(2253)/  -33.63500,  -32.88289,  -33.86051/
      DATA P(2254),P(2255),P(2256)/  -32.77917,  -33.41106,  -32.28661/
      DATA P(2257),P(2258),P(2259)/  -32.58559,  -31.15725,  -31.15697/
      DATA P(2260),P(2261),P(2262)/  -29.59799,  -29.11103,  -27.27611/
      DATA P(2263),P(2264),P(2265)/  -26.41639,  -24.41592,  -23.15949/
      DATA P(2266),P(2267),P(2268)/  -20.76700,  -18.99300,   -8.00000/
      DATA P(2269),P(2270),P(2271)/  -11.72448,  -12.72712,  -16.32320/
      DATA P(2272),P(2273),P(2274)/  -16.96900,  -20.19300,  -20.67400/
      DATA P(2275),P(2276),P(2277)/  -23.51900,  -23.69600,  -26.17800/
      DATA P(2278),P(2279),P(2280)/  -26.09800,  -28.44780,  -28.14500/
      DATA P(2281),P(2282),P(2283)/  -30.22300,  -29.68500,  -31.41300/
      DATA P(2284),P(2285),P(2286)/  -30.68051,  -32.06900,  -31.07075/
      DATA P(2287),P(2288),P(2289)/  -32.24661,  -31.07559,  -31.84326/
      DATA P(2290),P(2291),P(2292)/  -30.55735,  -30.97047,  -29.56330/
      DATA P(2293),P(2294),P(2295)/  -29.52023,  -27.67908,  -27.36207/
      DATA P(2296),P(2297),P(2298)/  -25.28345,  -24.70728,  -22.30359/
      DATA P(2299),P(2300),P(2301)/  -20.95985,  -16.22939,  -13.09700/
      DATA P(2302),P(2303),P(2304)/   -2.90500,   -4.44500,   -7.94900/
      DATA P(2305),P(2306),P(2307)/   -9.13500,  -12.19900,  -13.40400/
      DATA P(2308),P(2309),P(2310)/  -16.11800,  -16.99000,  -19.46800/
      DATA P(2311),P(2312),P(2313)/  -19.98200,  -22.19700,  -22.43000/
      DATA P(2314),P(2315),P(2316)/  -24.50800,  -24.41300,  -26.18800/
      DATA P(2317),P(2318),P(2319)/  -25.94900,  -27.43400,  -26.96400/
      DATA P(2320),P(2321),P(2322)/  -28.27100,  -27.46600,  -28.37684/
      DATA P(2323),P(2324),P(2325)/  -27.51047,  -28.11823,  -27.06419/
      DATA P(2326),P(2327),P(2328)/  -27.19611,  -25.99746,  -25.77528/
      DATA P(2329),P(2330),P(2331)/  -24.35983,  -23.83481,  -22.26706/
      DATA P(2332),P(2333),P(2334)/  -21.04440,  -16.76256,  -13.64720/
      DATA P(2335),P(2336),P(2337)/   -9.25386,   -3.06100,   -6.82217/
      DATA P(2338),P(2339),P(2340)/   -7.51700,  -10.99300,  -11.56900/
      DATA P(2341),P(2342),P(2343)/  -14.62300,  -14.87600,  -17.64200/
      DATA P(2344),P(2345),P(2346)/  -17.81000,  -20.32519,  -20.30700/
      DATA P(2347),P(2348),P(2349)/  -22.57700,  -22.28100,  -24.25000/
      DATA P(2350),P(2351),P(2352)/  -23.78000,  -25.42000,  -24.79600/
      DATA P(2353),P(2354),P(2355)/  -26.10000,  -25.23474,  -26.25363/
      DATA P(2356),P(2357),P(2358)/  -25.29324,  -25.94789,  -24.80057/
      DATA P(2359),P(2360),P(2361)/  -25.12355,  -23.78373,  -23.80060/
      DATA P(2362),P(2363),P(2364)/  -22.46707,  -21.76356,  -17.62871/
      DATA P(2365),P(2366),P(2367)/  -14.74263,  -10.49656,   -7.55675/
      DATA P(2368),P(2369),P(2370)/   -3.26000,    -.18802,   -2.13500/
      DATA P(2371),P(2372),P(2373)/   -3.27900,   -6.09400,   -7.29100/
      DATA P(2374),P(2375),P(2376)/   -9.77600,  -10.69500,  -12.99100/
      DATA P(2377),P(2378),P(2379)/  -13.62900,  -15.77900,  -16.07000/
      DATA P(2380),P(2381),P(2382)/  -17.93000,  -18.06400,  -19.62258/
      DATA P(2383),P(2384),P(2385)/  -19.53865,  -20.88693,  -20.36061/
      DATA P(2386),P(2387),P(2388)/  -21.45163,  -20.79606,  -21.54721/
      DATA P(2389),P(2390),P(2391)/  -20.67436,  -21.07534,  -20.04309/
      DATA P(2392),P(2393),P(2394)/  -20.06883,  -18.88446,  -18.27289/
      DATA P(2395),P(2396),P(2397)/  -14.80615,  -11.86897,   -8.13051/
      DATA P(2398),P(2399),P(2400)/   -5.23981,   -1.21219,    1.70682/
      DATA P(2401),P(2402),P(2403)/    5.77500,   -4.55500,   -4.98000/
      DATA P(2404),P(2405),P(2406)/   -7.89700,   -8.28900,  -10.91328/
      DATA P(2407),P(2408),P(2409)/  -11.13600,  -13.49500,  -13.44500/
      DATA P(2410),P(2411),P(2412)/  -15.51600,  -15.28100,  -17.01400/
      DATA P(2413),P(2414),P(2415)/  -16.57200,  -17.97500,  -17.31380/
      DATA P(2416),P(2417),P(2418)/  -18.34380,  -17.54432,  -18.19651/
      DATA P(2419),P(2420),P(2421)/  -17.15977,  -17.48315,  -16.37959/
      DATA P(2422),P(2423),P(2424)/  -15.96823,  -12.44768,  -10.38452/
      DATA P(2425),P(2426),P(2427)/   -6.66715,   -4.48426,    -.54529/
      DATA P(2428),P(2429),P(2430)/    1.77468,    5.82500,    8.35156/
      DATA P(2431),P(2432),P(2433)/     .17500,    -.95900,   -3.21000/
      DATA P(2434),P(2435),P(2436)/   -4.00300,   -6.25000,   -6.75100/
      DATA P(2437),P(2438),P(2439)/   -8.72500,   -9.04200,  -10.72437/
      DATA P(2440),P(2441),P(2442)/  -10.76003,  -12.25174,  -11.86578/
      DATA P(2443),P(2444),P(2445)/  -13.00705,  -12.48272,  -13.24970/
      DATA P(2446),P(2447),P(2448)/  -12.49831,  -12.89311,  -11.98711/
      DATA P(2449),P(2450),P(2451)/  -11.66155,   -8.63061,   -6.59391/
      DATA P(2452),P(2453),P(2454)/   -3.39398,   -1.26567,    2.24382/
      DATA P(2455),P(2456),P(2457)/    4.38698,    8.08673,   10.52260/
      DATA P(2458),P(2459),P(2460)/   14.25300,   16.89700,   20.80000/
      DATA P(2461),P(2462),P(2463)/   23.60400,    2.15100,    1.54700/
      DATA P(2464),P(2465),P(2466)/   -1.13644,   -1.57500,   -4.02700/
      DATA P(2467),P(2468),P(2469)/   -4.15900,   -6.31800,   -6.22600/
      DATA P(2470),P(2471),P(2472)/   -8.04400,   -7.76100,   -9.16600/
      DATA P(2473),P(2474),P(2475)/   -8.63791,   -9.65844,   -8.96411/
      DATA P(2476),P(2477),P(2478)/   -9.61314,   -8.76963,   -8.67323/
      DATA P(2479),P(2480),P(2481)/   -5.71159,   -4.33492,   -1.18375/
      DATA P(2482),P(2483),P(2484)/     .24047,    3.64638,    5.20362/
      DATA P(2485),P(2486),P(2487)/    8.82575,   10.60427,   14.39600/
      DATA P(2488),P(2489),P(2490)/   16.36679,   20.29700,   22.44000/
      DATA P(2491),P(2492),P(2493)/   26.49200,   28.77400,   32.98100/
      DATA P(2494),P(2495),P(2496)/   35.47500,    6.05400,    3.71300/
      DATA P(2497),P(2498),P(2499)/    3.06200,     .97600,     .55200/
      DATA P(2500),P(2501),P(2502)/   -1.24451,   -1.40946,   -2.92546/
      DATA P(2503),P(2504),P(2505)/   -2.66980,   -3.80476,   -3.35494/
      DATA P(2506),P(2507),P(2508)/   -4.16440,   -3.54435,   -3.56311/
      DATA P(2509),P(2510),P(2511)/    -.97365,     .30369,    2.96948/
      DATA P(2512),P(2513),P(2514)/    4.30020,    7.04519,    8.60779/
      DATA P(2515),P(2516),P(2517)/   11.46946,   13.26974,   16.34209/
      DATA P(2518),P(2519),P(2520)/   18.37904,   21.64374,   23.85268/
      DATA P(2521),P(2522),P(2523)/   27.33322,   29.65240,   33.27600/
      DATA P(2524),P(2525),P(2526)/   35.79300,   39.59800,   42.29600/
      DATA P(2527),P(2528),P(2529)/   46.25300,    8.57900,    6.03400/
      DATA P(2530),P(2531),P(2532)/    5.76300,    3.52300,    3.47200/
      DATA P(2533),P(2534),P(2535)/    1.65400,    1.81100,     .41600/
      DATA P(2536),P(2537),P(2538)/     .83275,    -.20168,     .32216/
      DATA P(2539),P(2540),P(2541)/     .08490,    2.51894,    3.27737/
      DATA P(2542),P(2543),P(2544)/    5.87401,    6.63591,    9.37901/
      DATA P(2545),P(2546),P(2547)/   10.26010,   12.95500,   14.30944/
      DATA P(2548),P(2549),P(2550)/   17.22997,   18.81804,   21.98741/
      DATA P(2551),P(2552),P(2553)/   23.66232,   27.17231,   28.93602/
      DATA P(2554),P(2555),P(2556)/   32.43491,   34.54424,   38.39600/
      DATA P(2557),P(2558),P(2559)/   40.70000,   44.70700,   47.08500/
      DATA P(2560),P(2561),P(2562)/   11.26500,   10.69800,    8.91322/
      DATA P(2563),P(2564),P(2565)/    8.62265,    7.12425,    7.27631/
      DATA P(2566),P(2567),P(2568)/    6.12335,    6.42086,    6.00891/
      DATA P(2569),P(2570),P(2571)/    8.12381,    8.69333,   10.82866/
      DATA P(2572),P(2573),P(2574)/   11.55511,   13.74150,   14.50871/
      DATA P(2575),P(2576),P(2577)/   16.60747,   17.81578,   20.22128/
      DATA P(2578),P(2579),P(2580)/   21.62983,   24.30253,   25.84614/
      DATA P(2581),P(2582),P(2583)/   28.89012,   30.67491,   33.55720/
      DATA P(2584),P(2585),P(2586)/   35.91049,   39.14368,   41.49800/
      DATA P(2587),P(2588),P(2589)/   45.10300,   47.60100,   51.39800/
      DATA P(2590),P(2591),P(2592)/   14.00100,   13.84000,   12.03200/
      DATA P(2593),P(2594),P(2595)/   12.07400,   10.66600,   10.92325/
      DATA P(2596),P(2597),P(2598)/   10.29391,   12.17106,   12.35882/
      DATA P(2599),P(2600),P(2601)/   14.45795,   14.65531,   16.92664/
      DATA P(2602),P(2603),P(2604)/   17.18991,   19.37056,   19.98916/
      DATA P(2605),P(2606),P(2607)/   22.30131,   23.18552,   25.80132/
      DATA P(2608),P(2609),P(2610)/   26.76307,   29.57991,   30.85720/
      DATA P(2611),P(2612),P(2613)/   33.81049,   35.44368,   38.72865/
      DATA P(2614),P(2615),P(2616)/   40.60894,   44.25008,   46.30500/
      DATA P(2617),P(2618),P(2619)/   50.20200,   52.39000,   19.73203/
      DATA P(2620),P(2621),P(2622)/   19.31847,   17.78931,   17.80052/
      DATA P(2623),P(2624),P(2625)/   17.03569,   18.63724,   18.51842/
      DATA P(2626),P(2627),P(2628)/   20.37782,   20.36594,   22.10000/
      DATA P(2629),P(2630),P(2631)/   22.32208,   23.86008,   24.32612/
      DATA P(2632),P(2633),P(2634)/   26.01919,   26.82065,   28.91072/
      DATA P(2635),P(2636),P(2637)/   29.89034,   32.16687,   33.42097/
      DATA P(2638),P(2639),P(2640)/   35.93864,   37.48353,   40.33585/
      DATA P(2641),P(2642),P(2643)/   42.32407,   45.34063,   47.63607/
      DATA P(2644),P(2645),P(2646)/   50.76367,   53.21600,   56.80200/
      DATA P(2647),P(2648),P(2649)/   21.87800,   23.20856,   23.01900/
      DATA P(2650),P(2651),P(2652)/   24.54600,   24.28400,   25.82376/
      DATA P(2653),P(2654),P(2655)/   25.70005,   27.37136,   27.32980/
      DATA P(2656),P(2657),P(2658)/   29.00678,   29.21757,   31.20145/
      DATA P(2659),P(2660),P(2661)/   31.60316,   33.80313,   34.60153/
      DATA P(2662),P(2663),P(2664)/   36.91343,   38.14058,   40.91407/
      DATA P(2665),P(2666),P(2667)/   42.44063,   45.38607,   47.30367/
      DATA P(2668),P(2669),P(2670)/   50.56873,   52.70927,   56.19700/
      DATA P(2671),P(2672),P(2673)/   58.61400,   31.57735,   32.72300/
      DATA P(2674),P(2675),P(2676)/   32.56341,   33.70100,   33.76374/
      DATA P(2677),P(2678),P(2679)/   35.22220,   35.61396,   37.35200/
      DATA P(2680),P(2681),P(2682)/   37.94194,   39.95043,   41.03778/
      DATA P(2683),P(2684),P(2685)/   43.37010,   44.86750,   47.45073/
      DATA P(2686),P(2687),P(2688)/   49.30527,   52.32092,   54.25604/
      DATA P(2689),P(2690),P(2691)/   57.41300,   59.87000,   63.20200/
      DATA P(2692),P(2693),P(2694)/   36.07460,   37.38918,   36.92964/
      DATA P(2695),P(2696),P(2697)/   38.43200,   38.35840,   40.04266/
      DATA P(2698),P(2699),P(2700)/   40.33804,   42.17944,   42.89351/
      DATA P(2701),P(2702),P(2703)/   45.08782,   46.15869,   48.58348/
      DATA P(2704),P(2705),P(2706)/   50.12132,   52.95104,   54.71301/
      DATA P(2707),P(2708),P(2709)/   57.74984,   59.79972,   63.09815/
      DATA P(2710),P(2711),P(2712)/   65.38941,   68.99600,   42.43900/
      DATA P(2713),P(2714),P(2715)/   43.39800,   43.28900,   44.52000/
      DATA P(2716),P(2717),P(2718)/   44.73900,   46.17400,   46.54745/
      DATA P(2719),P(2720),P(2721)/   48.41704,   49.38639,   51.50027/
      DATA P(2722),P(2723),P(2724)/   52.93023,   55.46398,   57.16829/
      DATA P(2725),P(2726),P(2727)/   59.87590,   61.89348,   64.98888/
      DATA P(2728),P(2729),P(2730)/   67.14800,   70.55600,   73.10400/
      DATA P(2731),P(2732),P(2733)/   47.32000,   46.79800,   48.05700/
      DATA P(2734),P(2735),P(2736)/   47.88300,   49.26800,   49.38436/
      DATA P(2737),P(2738),P(2739)/   51.18600,   51.71566,   53.69759/
      DATA P(2740),P(2741),P(2742)/   54.79916,   57.17719,   58.44784/
      DATA P(2743),P(2744),P(2745)/   60.99942,   62.61274,   65.52763/
      DATA P(2746),P(2747),P(2748)/   67.38636,   70.74422,   72.98319/
      DATA P(2749),P(2750),P(2751)/   76.64134,   79.05601,   52.70400/
      DATA P(2752),P(2753),P(2754)/   53.40300,   53.21400,   54.27400/
      DATA P(2755),P(2756),P(2757)/   54.36400,   55.65600,   56.09800/
      DATA P(2758),P(2759),P(2760)/   57.79900,   58.68558,   60.70348/
      DATA P(2761),P(2762),P(2763)/   61.80964,   63.96274,   65.48266/
      DATA P(2764),P(2765),P(2766)/   68.07401,   69.84335,   72.94578/
      DATA P(2767),P(2768),P(2769)/   75.22134,   78.52801,   80.92900/
      DATA P(2770),P(2771),P(2772)/   84.39301,   57.81800,   57.20300/
      DATA P(2773),P(2774),P(2775)/   58.29200,   58.02700,   59.35100/
      DATA P(2776),P(2777),P(2778)/   59.32565,   60.93900,   61.46965/
      DATA P(2779),P(2780),P(2781)/   63.37800,   64.08567,   66.12865/
      DATA P(2782),P(2783),P(2784)/   67.23344,   69.71935,   71.16609/
      DATA P(2785),P(2786),P(2787)/   74.12834,   76.02814,   79.29509/
      DATA P(2788),P(2789),P(2790)/   81.33450,   84.80300,   87.03900/
      DATA P(2791),P(2792),P(2793)/   64.19901,   63.95900,   64.92400/
      DATA P(2794),P(2795),P(2796)/   64.86100,   66.10700,   66.43201/
      DATA P(2797),P(2798),P(2799)/   67.96600,   68.60400,   70.28900/
      DATA P(2800),P(2801),P(2802)/   71.17001,   73.26601,   74.50423/
      DATA P(2803),P(2804),P(2805)/   77.28815,   79.00742,   81.98640/
      DATA P(2806),P(2807),P(2808)/   84.08259,   87.18001,   89.40301/
      DATA P(2809),P(2810),P(2811)/   68.40000,   69.40601,   69.00201/
      DATA P(2812),P(2813),P(2814)/   70.21201,   70.12439,   71.55600/
      DATA P(2815),P(2816),P(2817)/   71.89680,   73.61100,   74.06751/
      DATA P(2818),P(2819),P(2820)/   75.97867,   76.81105,   79.34117/
      DATA P(2821),P(2822),P(2823)/   80.89819,   83.79296,   85.47996/
      DATA P(2824),P(2825),P(2826)/   88.58380,   90.41901,   93.69701/
      DATA P(2827),P(2828),P(2829)/   75.46700,   76.32001,   76.20000/
      DATA P(2830),P(2831),P(2832)/   77.22900,   77.31600,   78.70000/
      DATA P(2833),P(2834),P(2835)/   79.10101,   80.69501,   81.30100/
      DATA P(2836),P(2837),P(2838)/   83.57800,   84.83599,   87.60957/
      DATA P(2839),P(2840),P(2841)/   88.98994,   91.68258,   93.61700/
      DATA P(2842),P(2843),P(2844)/   96.54501,   81.80701,   81.49900/
      DATA P(2845),P(2846),P(2847)/   82.86600,   82.87120,   84.43900/
      DATA P(2848),P(2849),P(2850)/   84.71820,   86.84546,   87.81741/
      DATA P(2851),P(2852),P(2853)/   90.21777,   91.47301,   94.10300/
      DATA P(2854),P(2855),P(2856)/   95.60500,   98.49901,  100.15400/
      DATA P(2857),P(2858),P(2859)/   87.89600,   88.79900,   88.73200/
      DATA P(2860),P(2861),P(2862)/   89.97100,   90.14001,   91.99700/
      DATA P(2863),P(2864),P(2865)/   92.78201,   94.90301,   95.93501/
      DATA P(2866),P(2867),P(2868)/   98.34000,   99.61501,  102.17700/
      DATA P(2869),P(2870),P(2871)/  103.76200,   93.78201,   93.30400/
      DATA P(2872),P(2873),P(2874)/   94.53900,   94.24815,   96.01100/
      DATA P(2875),P(2876),P(2877)/   96.47301,   98.39201,   99.14201/
      DATA P(2878),P(2879),P(2880)/  101.30200,  102.38800,  104.83000/
      DATA P(2881),P(2882),P(2883)/  106.17000,  100.04100,  100.70400/
      DATA P(2884),P(2885),P(2886)/  100.46900,  101.94100,  102.20500/
      DATA P(2887),P(2888),P(2889)/  103.79400,  104.42600,  106.33300/
      DATA P(2890),P(2891),P(2892)/  107.19400,  109.42500,  110.53000/
      DATA P(2893),P(2894),P(2895)/  105.39900,  106.79800,  106.59590/
      DATA P(2896),P(2897),P(2898)/  108.23800,  108.49800,  110.20800/
      DATA P(2899),P(2900),P(2901)/  110.77700,  112.77200,  113.57500/
      DATA P(2902),P(2903),P(2904)/  113.45900,  113.45600,  114.58200/
      DATA P(2905),P(2906),P(2907)/  114.71000,  116.18600,  116.62100/
      DATA P(2908),P(2909),P(2910)/  118.30800,  118.98900,  119.89300/
      DATA P(2911),P(2912),P(2913)/  119.61150,  121.09600,  121.13300/
      DATA P(2914),P(2915),P(2916)/  122.74600,  123.10200,  124.92700/
      DATA P(2917),P(2918),P(2919)/  127.21100,  128.49000,  128.10500/
      DATA P(2920),P(2921),P(2922)/  129.30600,  129.57600,  131.08300/
      DATA P(2923),P(2924),P(2925)/  131.55400,  134.09400,  133.69600/
      DATA P(2926),P(2927),P(2928)/  135.20000,  134.71800,  136.07000/
      DATA P(2929),P(2930),P(2931)/  136.28700,  139.02100,  142.96300/
*
      DATA IZ/   2,   8,  16,  25,  35,  48,  63,  78,  93, 109, 126,
     * 144, 162, 181, 202, 225, 249, 273, 297, 321, 345, 369, 393, 417,
     * 441, 465, 490, 515, 544, 573, 602, 631, 660, 690, 718, 746, 775,
     * 807, 839, 869, 899, 929, 960, 991,1023,1056,1089,1123,1158,1195,
     *1233,1270,1307,1344,1382,1422,1462,1501,1540,1579,1615,1651,1687,
     *1723,1757,1791,1825,1859,1893,1927,1961,1996,2029,2062,2095,2128,
     *2163,2198,2233,2268,2302,2336,2370,2402,2431,2462,2495,2528,2560,
     *2590,2619,2647,2672,2692,2712,2731,2751,2771,2791,2809,2827,2843,
     *2857,2870,2882,2893,2902,2910,2917,2924,2931,2932/
*
      DATA IA/   1,   3,   4,   5,   7,   8,  10,  12,  14,  16,  18,
     *  20,  21,  22,  24,  26,  28,  30,  32,  34,  36,  38,  40,  42,
     *  44,  45,  48,  50,  52,  54,  56,  58,  60,  65,  67,  69,  71,
     *  73,  77,  79,  81,  83,  85,  87,  89,  91,  94,  96,  98, 100,
     * 103, 106, 108, 110, 112, 114, 117, 119, 121, 126, 128, 130, 132,
     * 136, 138, 140, 142, 144, 146, 148, 150, 154, 156, 158, 160, 162,
     * 165, 168, 171, 175, 177, 181, 185, 190, 193, 196, 200, 203, 207,
     * 210, 213, 218, 225, 228, 231, 233, 235, 237, 240, 242, 245, 249,
     * 251, 253, 255, 258, 260, 263, 265, 267, 272,   0/
*
      END
************************************************************************
* Readme.txt
*===========
*
* Copy the rest of this file and paste it to an input for this code. 
* The lines starting with a `*' are commented and the code
* skips reading them. Include as many commented lines as desired
* to remember how to construct the input and keep track of changes in
* the input.  
* 
* To start with, delete the `*' from the three five input data rows, 
* and calculate the bound-state wavefunction. As a next step, delete
* the '*' from the following input rows, keeping the `*' whenever the
* input is not used (depends upon the options IPOT, IDPROJ and IDTARG).
*
*    ********    Input of subroutine MOMDIS   *******
*
* Read write-up for more details.
*=========================================
* General input
*==============
*
* ZPROJ, APROJ = projectile charge and mass numbers
* ZTARG, ATARG = target charge and mass numbers
* ECA = laboratory energy per nucleon (MeV)
*
* ZPROJ       APROJ       ZTARG       ATARG      ECA
*  6           15           4           9         103.
*
* ZCORE,  ACORE  = projectile core charge and mass number
*
* ZCORE       ACORE
*  6           14
*
* S2C = Spectroscopic factor
*
*  S2C
*   1
*
*=========================================
* Input for bound state function
*===============================
*
* If your option is to calculate the bound-state wavefunction,
* enter the next two input lines, otherwise comment them with a'*' 
*
* N_0  =  nodes of the Wave Function (exclude origin)   
* J0   =  total angular momentum   
* L0   =  orbital angular momentum
*
* N_0    J0     L0    
* 
*  1     0.5     0     
*
* V0  =  depth of central pot. 
* VS0 =  depth of spin-orbit potential 
* R0  =  radius of the potential 
* AA  =  diffuseness of the pot. 
* RS0 =  radius of the spin-orbit pot. 
* AAS =  diffuseness of the spin-orbit pot. 
* RC  =  Coulomb radius (usually, RC = R0)  
*-----------------------------------------------------------------
*  WS =  V_0 f(r,R0,AA) - V_S0 (l.s) (r_0^2/r) d/dr f(r,RS0,AAS) 
*  
*  f(r,R0,a) = [ 1 + exp((r-R_0)/a) ]^(-1)
*
*  r_0 = 1.4138 fm is the pion Compton wavelength. 
*-----------------------------------------------------------------
* V0        R0       AA     VS0     RS0      AAS     RC
*
* -61.849     2.67      0.6    0.0    2.4    0.6       2.67 
*
*=========================================
* Input for S-matrices
*=====================
*
* Comment the two previous input data lines with a'*' 
*
* ISMAT = 0, or 1
* If ISMAT=0, then S-matrix is for nucleon+target
* If ISMAT=1, then S-matrix is for core+target
* ISMAT
*  0
*  1
*
* IPOT=0, 1, or 2 
* If IPOT=0, then optical potential is a Woods-Saxon, with real and 
* imaginary parts 
* If IPOT=1, then optical potential is entered in a file
* If IPOT=2, then optical potential is built with the 't-rho-rho' method
* IPOT
*  2   
*
* If IPOT=0, enter V0_ws [VI_ws]  = real part [imaginary] (>0, both) of Woods-Saxon.
*                  r0_ws [r0I_ws] = radius parameter (R_ws = r0 * (ap^1/3 + at^1/3)).
*                  d_ws [dI_ws]   = diffuseness.
* If IOPW is not equal to 0, place a '*' sign at the beginning of this line, or delete it.
*
*   V0 [MeV]     r0[fm]      d[fm]      VI [MeV]     r0_I [fm]     dI [fm]
*    50.          1.067       0.8        58.          1.067         0.8
*
* If IPOT=1, enter name of fileS where optical potential is located.
* In the input files the first row must give the number of the remaining rows 
* which form an ordered list of  R (in fm) x V in steps of a constant value of DR. 
*
*  FILENAME1          FILENAME2
*  realpot.inp        imagpot.inp
* 
* If IPOT=2, enter the following options
* IPAULI = 0 for no Pauli correction of nn cross sections
*        = 1 for Pauli correction of nn cross sections
*
* IPAULI
* 0
*
* Enter options for projectile and target densities  
* IDPROJ,IDTARG=1,2,3,4,5 
* If IDPROJ=1, projectile density is a Gaussian times (1+beta*(r/alpha)^pwr)
* If IDPROJ=2, projectile density is a Yukawa times r^pwr
* If IDPROJ=3, projectile density is a Fermi function times (1+w*(r/c)^pwr)
* If IDPROJ=4, projectile density is a (Fermi function)^pwr
* If IDPROJ=5, projectile density is a calculated from liquid-drop model
* If IDPROJ=10, projectile density is from an input file
* The same applies for the target density
*
* IDPROJ      IDTARG   
* 1           1  
* 
* If IDPROJ=1 or 2, enter Gaussian or Yukawa density parameters for the projectile
*    =============
* Gaussian: (1+beta*(r/alpha)^pwr) * exp(-(r/alpha)**2)
* Yukawa: r**pwr * exp(-alpha/r)
*
* ALPHA=range, BETA=power multiplier, PWR=power
* ALPHA        BETA     PWR
*   0.7         0.       0.
*  1.73        1.38      1.
* 
* If IDPROJ=3, enter Fermi density parameters for the projectile
*    ======== 
* [1 + w * (r/ar)**pwr] / [1 + exp((r-rr)/ar)]
*
* RR=range,  AR=diffuseness,  W=coefficient of r^pwr,  PWR=power 
* RR          AR         W          PWR
* 2.          1.         1.         1.
* 
* If IDPROJ=4, enter Fermi^pwr density parameters for the projectile
*    ========
* 1/[1 + exp(r-rr)/ar)]**pwr 
* RR=range, AR=diffuseness, PWR=power, IDER=0(1): if IDER=0 use derivative of Fermi^pwr 
* RR         AR        PWR      IDER
* 2.         1.        1.       0
*
* If IDPROJ=10, enter  name of projectile file density 
*    =========
* ALPHAP is the proton size correction 
* FILENAME      ALPHAP
*  projd.inp     0.6
* 
* If IDTARG=1 or 2, enter Gaussian or Yukawa density parameters for the target
*    =============
* ALPHA=range, PWR=power
* ALPHA         BETA        PWR
* 1.93          0.          0.
* 
* If IDTARG=3, enter Fermi density parameters for the target
*    ========
* RR=range,  AR=diffuseness,  W=coefficient of r^pwr,  PWR=power 
* RR          AR         W          PWR
* 2.          1.         1.         1.
* 
* If IDTARG=4, enter Fermi^pwr density parameters for the target
*    ========
* RR=range, AR=diffuseness, PWR=power, IDER=0(1): if IDER=0 use derivative of Fermi^pwr 
* RR         AR        PWR      IDER
* 2.         1.        1.       0
*
* If IDTARG=10, enter  name of target file density 
*    =========
* ALPHAP is the proton size correction.
* In the input files the first row must give the number of the remaining rows 
* which form an ordered list of  R (in fm) x RHO (in fm^-3) in steps of a constant 
* value of DR. 
*
* FILENAME      ALPHAP
*  targd.inp     0.6
*
*==================
* END of Readme.txt
*==================
*********************************************************************