* The same as VLADCS but you can vary Rws and Aws for s1/2,p1/2,d5/2 * * IMPLICIT REAL*8(A-H,O-Z) REAL*4 RGRIE,RGRICS CHARACTER infile*100,outfile*100 COMPLEX*16 SLPL,SLMI,SCOUL,RRJ COMMON * /SAMS/PCENTR(250,6),PSL(250,2) * /SVO/SD(15,2) DIMENSION PHPL(501),PHMI(501),PHCOUL(501),PL(501),PL1(501) DIMENSION CSC(10000),poten(3,2) DATA DMH/0.04823D0/,DME/.06939D0/ WRITE(*,*) ' Name of input-file: ' READ(*,'(A)') infile WRITE(*,*) ' Name of output-file: ' READ(*,'(A)') outfile OPEN( 8,FILE=outfile,ACCESS='sequential',STATUS='unknown') OPEN(10,FILE='arrays.out',ACCESS='sequential',STATUS='unknown') OPEN( 9,FILE=infile,ACCESS='sequential',STATUS='unknown') RRJ=CMPLX(0.D0,1.D0) READ(9,*) AP READ(9,*) ZP READ(9,*) AT READ(9,*) ZT READ(9,*) PRM IF(DABS(PRM).LT.1D-7) PRM=AP*AT/(AP+AT) C *** for exact calculations factor of reduced mass should be given in C *** inputcs.dat file as M_p*M_t/(M_p+M_t), where all masses are in MeV, C *** divided by 938.92 MeV (i.e. by the average nucleon mass used in the C *** code corresponding to the factor DMH = 0.04823) READ(9,*) UWS READ(9,*) ROWS READ(9,*) DIFWS READ(9,*) USO READ(9,*) ROSO READ(9,*) DIFSO READ(9,*) ROCUL READ(9,*) THETA C *** theta = cm scattering angle (in degrees) READ(9,*) LMIN READ(9,*) LMAX C *** only partial waves from lmin to lmax are to be included READ(9,*) EMIN READ(9,*) EMAX READ(9,*) HENERG C *** 'Emin' and 'Emax' give an interval of the cm energy in which C *** the calculations are to be done with step of 'henerg'. READ(9,*) KEYC C *** if keyc not equal 1 then Coulomb amplitude is neglegted (not phases) READ(9,*) IPRINT C *** if iprint not equal 1 then no print potentials, polinomials, etc. C*** SD(I,K): I=s1/2,p1/2,p3/2,... K=n,p (potential depth factors) DO 8 I=1,15 SD(I,1)=1.D0 SD(I,2)=1.D0 8 CONTINUE C *** initially, all SD-factors = 1 !!! READ(9,*) NSD DO 9 I=1,NSD READ(9,*) SD(I,2) C write(8,16) i,SD(i,2) 16 FORMAT('# pot. depth factor',I2, ' =',F7.4) 9 CONTINUE DO I=1,3 READ (9,*) poten(i,1),poten(i,2) ENDDO ROWSs=poten(1,1) DIFWSs=poten(1,2) ROWSp=poten(2,1) DIFWSp=poten(2,2) ROWSd=poten(3,1) DIFWSd=poten(3,2) PRINT *,ROWSs,DIFWSs PRINT *,ROWSp,DIFWSp PRINT *,ROWSd,DIFWSd C *** the potential depth factors for s1/2, p1/2, p3/2, d3/2, d3/2, ... C *** total number of SD-factors = nsd C write(8,15) lmin,lmax,Emin,Emax,henerg 15 FORMAT('# lmin, lmax = =',2I4/, *'# Emin, Emax, he = = =',3F10.3,' (all in MeV)') C write(8,14) AP,ZP,AT,ZT,PRM,theta 14 FORMAT('# A_p, Z_p = =',2F9.1/, *'# A_t, Z_t = =',2F9.1/, *'# M_red =',F10.4/, *'# scattering angle =',F8.2,' degrees'/, *'# c.m. energy (MeV), cross section (mb/sr)') C _____ Preparation of the potential _______ RWS=ROWS*AT**(1.D0/3.D0) RWSs=ROWSs*AT**(1.D0/3.D0) RWSp=ROWSp*AT**(1.D0/3.D0) RWSd=ROWSd*AT**(1.D0/3.D0) RSO=ROSO*AT**(1.D0/3.D0) RCU=ROCUL*AT**(1.D0/3.D0) NPOT=250 HR=0.1D0 PILAM2=2.D0 C pion Compton wavelength squared pilam2=(hbar/m_pi*c)**2 = 2 fm**2 USOC=-USO*PILAM2/DIFSO ESQ=DME/DMH ZAV=ZP*ZT VCUL1=0.5D0*ZAV*ESQ/RCU VCUL2=ZAV*ESQ DO 91 K=1,NPOT R=HR*K PCENTR(K,2)=UWS/(1.D0+DEXP((R-RWS)/DIFWS)) PCENTR(K,4)=UWS/(1.D0+DEXP((R-RWSs)/DIFWSs)) PCENTR(K,5)=UWS/(1.D0+DEXP((R-RWSp)/DIFWSp)) PCENTR(K,6)=UWS/(1.D0+DEXP((R-RWSd)/DIFWSd)) EXPSO=DEXP((R-RSO)/DIFSO) PSL(K,2)=(USOC/R)*EXPSO/(1.D0+EXPSO)**2 IF(R.LE.RCU) PCENTR(K,3)=VCUL1*(3.D0-(R/RCU)**2) IF(R.GT.RCU) PCENTR(K,3)=VCUL2/R 91 CONTINUE C------------------------print potentials--------- IF(IPRINT.NE.1) GO TO 103 WRITE(10,102) 102 FORMAT('Potentials: r(fm), pot_centr, pot_coul, pot_so (MeV)') DO 101 K=1,100 R=HR*K WRITE(10,100) R,PCENTR(K,2),PCENTR(K,3),PSL(K,2) 100 FORMAT(2X,F7.2,1P,3D15.7) 101 CONTINUE 103 CONTINUE C------------------------------end print--------------- LI=LMIN+1 LF=LMAX+1 DO 27 K=1,LF PHPL(K)=0.D0 PHMI(K)=0.D0 PHCOUL(K)=0.D0 PL(K)=0.D0 PL1(K)=0.D0 27 CONTINUE C *** calculations of the Legandre polinomials P_l and associated C *** Legandre polinomials P_l^1 for l from 0 to lmax PI=ACOS(-1.D0) X=DCOS(THETA*PI/180.D0) PLX=PLEG(X,PL,LMAX,0) PLY=PLEG(X,PL1,LMAX,1) C------------------print Legandre polinomials--------- IF(IPRINT.NE.1) GO TO 31 WRITE(10,21) X 21 FORMAT('Legandre polinomials P_l and P_l^1 at x =',F6.1) LF=LMAX+1 DO 22 KK=1,LF L=KK-1 WRITE(10,23) L,PL(KK),PL1(KK) 23 FORMAT('l =',I3,' P_l =',F10.5,' p_l^1 =',F10.5) 22 CONTINUE 31 CONTINUE C-----------------------------end print-------------- IT=2 REDMAS=PRM NE=(EMAX-EMIN)/HENERG+1 C *** calculations of Coulomb and nuclear phases C *** at given E for j=l+1/2 and j=l-1/2 for all l from lmin to lmax DO 4 K=1,NE E=EMIN+HENERG*(K-1) RK=DSQRT(DMH*PRM*E) ETA=0.5D0*DME*PRM*ZAV/RK PHC=SIGMAL(ETA,LMAX,PHCOUL) DO 5 LT=LI,LF RL=DFLOAT(LT-1) RJ=RL+0.5D0 WRITE(10,26) L,PHPL(KK),PHMI(KK),PHCOUL(KK) CALL PREPE(IT,ZAV,RJ,RL,E,PHASE,REDMAS) PHPL(LT)=PHASE RJ=RL-0.5D0 IF(RJ.LT.0.D0) GO TO 6 CALL PREPE(IT,ZAV,RJ,RL,E,PHASE,REDMAS) PHMI(LT)=PHASE 6 CONTINUE 5 CONTINUE C-----------print phases------------------- IF(IPRINT.NE.1) GO TO 32 WRITE(10,24) E 24 FORMAT('phases at E_cm =',F9.5) DO 25 KK=1,LF L=KK-1 WRITE(10,26) L,PHPL(KK),PHMI(KK),PHCOUL(KK) 26 FORMAT('l =',I3,' phase^+ =',F10.5,' phase^- =',F10.5, *' phase_coul =',F10.5) 25 CONTINUE 32 CONTINUE WRITE(10,34) 34 FORMAT(' Elements of S-matrix') DO 33 KK=1,LF LT=KK-1 SLPL= EXP(2.D0*RRJ*PHPL(KK)*PI/180.D0) SLMI= EXP(2.D0*RRJ*PHMI(KK)*PI/180.D0) SCOUL=EXP(2.D0*RRJ*PHCOUL(KK)) WRITE(10,35) LT,SLPL,SLMI,SCOUL 35 FORMAT('l =',I3,' S_l^+,S_l^-,S_l^Coul==='/,6D12.4) 33 CONTINUE WRITE(10,36) E,RK,ETA 36 FORMAT(' E, k, eta ===',3D15.7) C----------------end print----------- CSC(K)=DSDOM(THETA,ETA,RK,LMIN,LMAX,PHPL,PHMI,PHCOUL, *KEYC,PL,PL1) RGRIE=SNGL(E) RGRICS=SNGL(CSC(K)) WRITE(8,13) RGRIE,RGRICS 13 FORMAT(8X,F12.4,F15.7) 4 CONTINUE STOP END c-------------------------------------------------------------- SUBROUTINE PREPE(IT,ZAV,RJE,RLE,EE,PHASE,REDMAS) c-------------------------------------------------------------- IMPLICIT REAL*8(A-H,O-Z) C Calculates scattering wave function C with RJE,RLE,EE quantum numbers in the mean field C given by PT(250,k=2,3) and SL(250,k=2) (IT=2 - fixed) C /FSCAT/FSCT(5000) (from 0.1 to 500 fm) C normalized to delta(e-e') C*** BEFORE CALLING THIS SUBROUTINE THE POTENTIALS SHOULD BE C*** IN THE COMMON BLOCK /SAMS/ COMMON * /SAMS/PT(250,6),SL(250,2) * /PARFP/AA,ZZ,T,V1,V2,ROO,ROC,DIF,RLAM,RJ,RL * /FORFP/RN,RC,VO,VCO,VCD,RL1,XL,EX * /FSCAT/FSCT(5000) * /PRMAS/PRM AA=777.D0 ZZ=ZAV PRM=REDMAS RJ=RJE RL=RLE T=IT*1.D0-1.5D0 CALL PRPOT CALL WFSCAT(EE,PHASE) RETURN END C________________________________________________________ DOUBLE PRECISION FUNCTION SIGMAL(ETA,L,SLC) C *** sigmal(eta,l)= Coulomb phase = arg[Gamma(l+1+i*eta)] C *** all Coulomb phase from l=0 to given l are placed in slc(l) C *** so that sigma_l=slc(l+1) (in radians) IMPLICIT REAL*8(A-H,O-Z) DIMENSION SLC(2) DATA C/0.577215664901533D0/ DATA EPS/1.D-13/ SIGMA0=-C*ETA RL=DFLOAT(L) RM=0.D0 2 CONTINUE C1M=ETA/(1.D0+RM) CM=C1M-DATAN(C1M) IF(DABS(CM).LT.EPS) GO TO 1 SIGMA0=SIGMA0+CM RM=RM+1.D0 GO TO 2 1 CONTINUE SLC(1)=SIGMA0 SIGMAL=SIGMA0 IF(L.EQ.0) RETURN DO 3 M=1,L SLC(M+1)=SLC(M)+DATAN(ETA/DFLOAT(M)) 3 CONTINUE SIGMAL=SLC(L+1) RETURN END C________________________________________________ DOUBLE PRECISION FUNCTION PLEG(X,PLM,L,M) C *** associated Legendre polinomials P_l^m(x), L should be > or = M, C *** for l from M to L at fixed M, x=cos(Theta) C *** the polinomial P_n^m is placed in plm(n+1), where n<=L, C *** so that P_L^M = plm(L+1) IMPLICIT REAL*8(A-H,O-Z) DIMENSION PLM(2) IF(L.EQ.0.AND.M.EQ.0) THEN PLEG=1.D0 PLM(1)=1.D0 RETURN ENDIF IF(M.EQ.0) THEN PLM(1)=1.D0 PLM(2)=X PLEG=X IF(L.EQ.1) RETURN KMAX=L+1 DO 1 K=3,KMAX RLT=DFLOAT(K-1) PLM(K)=((2.D0*RLT-1.D0)*X*PLM(K-1)-(RLT-1.D0)*PLM(K-2))/RLT 1 CONTINUE PLEG=PLM(L+1) RETURN ENDIF ST=0.D0 IF(X.LT.1.D0) ST=DSQRT(1.D0-X**2) DF=1.D0 DO 2 N=1,M TM=DFLOAT(2*N-1) DF=DF*TM 2 CONTINUE DO 3 K=1,M PLM(K)=0.D0 3 CONTINUE PLM(M+1)=DF*ST**M IF(L.EQ.M) THEN PLEG=PLM(L+1) RETURN ENDIF KMIN=M+2 KMAX=L+1 RM=DFLOAT(M) DO 4 K=KMIN,KMAX RLT=DFLOAT(K-1) PLM(K)=((2.D0*RLT-1.D0)*X*PLM(K-1)-(RLT-1.D0+RM)*PLM(K-2)) */(RLT-RM) 4 CONTINUE PLEG=PLM(L+1) RETURN END C__________________________________________________________ DOUBLE PRECISION FUNCTION DSDOM(THETA,ETA,RK, *LMIN,LMAX,PHPL,PHMI,PHCOUL,KEYC,PL,PL1) C *** differential cross section dsigma/domega (in mb/sr) C *** theta - scattering angle (in degrees) C *** eta - Coulomb parameter (= Z*e^2/hbar*v = Z*e^2*k/2E = Ze^2*M/hbar^2*k) C *** rk - relative moution momentum [= sqrt(2ME/hbar^2)] C *** lmin - minimum angular momentum C *** lmax - maximum angular momentum C *** phpl - set of nuclear phases with j=l+1/2 C *** phmi - set of nuclear phases with j=l-1/2 C *** phcoul - set of Coulomb phases C *** keyc - if 0 - no Coulomb amplitude, if 1 - yes C *** pl - set of Legandre polinomials P_l C *** pl1 - set of associated Legandre polinomials P_l^1 IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 RJ,A,B,SLPL,SLMI,SCOUL,FCOUL,ARG,DIK DIMENSION PHPL(2),PHMI(2),PHCOUL(2),PL(2),PL1(2) RJ=CMPLX(0.D0,1.D0) DIK=-0.5D0*RJ/RK LI=LMIN+1 LF=LMAX+1 PI=ACOS(-1.D0) X=THETA*PI/180.D0 ST2=SIN(0.5D0*X) ARG=-2.D0*RJ*(ETA*LOG(ST2)-PHCOUL(1)) FCOUL=-0.5D0*ETA/RK/ST2**2*EXP(ARG) A=CMPLX(0.D0,0.D0) B=CMPLX(0.D0,0.D0) DO 1 LT=LI,LF RL=DFLOAT(LT-1) SLPL= EXP(2.D0*RJ*PHPL(LT)*PI/180.D0) SLMI= EXP(2.D0*RJ*PHMI(LT)*PI/180.D0) SCOUL=EXP(2.D0*RJ*PHCOUL(LT)) A=A+((RL+1.D0)*(SLPL-1.D0)+RL*(SLMI-1.D0))*PL(LT)*SCOUL B=B+(SLPL-SLMI)*PL1(LT)*SCOUL 1 CONTINUE A=DIK*A B=DIK*B IF(KEYC.EQ.1) A=A+FCOUL DSDOM=(ABS(A)**2+ABS(B)**2)*10.D0 RETURN END SUBROUTINE PRPOT IMPLICIT REAL*8(A-H,O-Z) COMMON /PARFP/A,Z,T,V1,V2,RO,ROC,DIF,RLAM,RJ,RL * /FORFP/RN,RC,VO,VCO,VCD,RL1,XL,EX * /FGP/FP1(250) * /PRMAS/PRM * /SVO/SD(15,2) * /SAMS/PT2(250,6),S2(250,2) DATA DMH/0.04823D0/,DME/.06939D0/,H/.1D0/,N/250/ IND=RL+RJ+1.D0 IT =T+1.6D0 IT1=T+1.6D0 RL1=RL*(RL+1.0D0) DMHPRM=DMH*PRM SP =DMHPRM IF(IND.LE.15) SP=SD(IND,IT)*DMHPRM XL =(RJ*(RJ+1.D0)-RL1-.75D0)*DMHPRM IF ((RL.EQ.0.).AND.(RJ.EQ.0.5)) IT1=4 IF ((RL.EQ.1.).AND.(RJ.EQ.0.5)) IT1=5 IF ((RL.EQ.2.).AND.(RJ.EQ.2.5)) IT1=6 DO 7 K=1,N R=K*H FP1(K)=SP*PT2(K,IT1)+XL*S2(K,IT)+RL1/(R*R) 7 CONTINUE VCO=0.D0 IF(IT.EQ.1) GO TO 4 VCO=DME*PRM*Z C *** here Z=ZAV=ZP*ZT when calling this 'PRPOT' subroutine DO 17 K=1,N FP1(K)=FP1(K)+DMHPRM*PT2(K,3) 17 CONTINUE 4 CONTINUE RETURN END DOUBLE PRECISION FUNCTION FPOTS(R) IMPLICIT REAL*8(A-H,O-Z) COMMON /FORFP/RN,RC,VO,VCO,VCD,RL1,XL,EX FPOTS=RL1/(R*R)+VCO/R RETURN END DOUBLE PRECISION FUNCTION BESR(RL,W) IMPLICIT REAL*8(A-H,O-Z) DATA D / 1.D-7 / BESR = 1.D0 B = 1.D0 AK = 0.D0 1 AK = AK+1.D0 B = -B*W/(AK*(AK+RL+0.5D0)) BESR = BESR+B IF (DABS(B).GT.D) GO TO 1 RETURN END c--------------------------------------------------------------- SUBROUTINE WFSCAT(E,PHASE) C********************************************************************** C Calculate scattering radial wave function, normalized C to delta(e-e'), and phase shift PHASE (in degrees) C the wave fuction is /FSCAT/FUN(5000) with step 0.1 fm C********************************************************************** IMPLICIT REAL*8(A-H,O-Z) DIMENSION F(51),FPR(51),G(51),GPR(51) COMMON * /PARFP/A,Z,T,V1,V2,RO,ROC,DIF,RLAM,RJ,RL * /FGP/FP1(250) * /PRMAS/PRM * /FSCAT/FUN(5000) DATA RMAX1/500.D0/,ACCUR/1.D-7/,STEP/100.D0/,NMAX/5000/ PI=ACOS(-1.D0) H=0.1D0 B11=H*H/12.D0 B22=B11*10.D0 P1=0.04823D0*PRM P2=0.06939D0*PRM E1=P1*E L1=RL+1.01D0 L=L1 IF(L1.GT.3) L=L1/2 RK2=0.25D0*(E1-FP1(L)+RL*(RL+1.D0)/((L*H)**2)) L2=L1+1 L3=L2+1 C1=1.D-4/(10.D0**L1) DO 1 K=1,L2 R=K*H FUN(K)=C1*R**L1*BESR(RL,RK2*R*R) 1 CONTINUE F1=FP1(L1)-E1 F2=FP1(L2)-E1 DO 3 K=L3,250 F3=FP1(K)-E1 FUN(K)=((2.D0+B22*F2)*FUN(K-1)- *(1.D0-B11*F1)*FUN(K-2))/(1.D0-B11*F3) F1=F2 F2=F3 3 CONTINUE C C FUN IS OBTAINED FROM 0.1 TO 25.0 FM (STEP H=0.1 FM) C DO 178 K=251,NMAX F3=FPOTS(K*H)-E1 FUN(K)=((2.D0+B22*F2)*FUN(K-1)- *(1.D0-B11*F1)*FUN(K-2))/(1.D0-B11*F3) F1=F2 F2=F3 178 CONTINUE C C FUN IS OBTAINED FROM 25.0 TO 500 FM (STEP H=0.1 FM) C NM1=NMAX-1 NM2=NMAX Y1S=FUN(NM1) Y2S=FUN(NM2) RK=DSQRT(E1) ETA=(T+0.5D0)*P2*0.5D0*Z/RK RTURN=(ETA+DSQRT(ETA*ETA+RL*(RL+1.D0)))/RK RM=RTURN+1.D0 IF(RM.GT.RMAX1) THEN NM=RM/H+3.01D0 NMAX1=NMAX+1 DO 77 K=NMAX1,NM F3=FPOTS(K*H)-E1 Y3S=((2.D0+B22*F2)*Y2S- *(1.D0-B11*F1)*Y1S)/(1.D0-B11*F3) F1=F2 F2=F3 Y1S=Y2S Y2S=Y3S 77 CONTINUE NM1=NM-1 NM2=NM ENDIF L=RL+0.1D0 ROG1=NM1*H*RK CALL RCWFN(ROG1,ETA,L,L,F,FPR,G,GPR,ACCUR,STEP) FC1=F(L+1) GC1=G(L+1) ROG2=NM2*H*RK CALL RCWFN(ROG2,ETA,L,L,F,FPR,G,GPR,ACCUR,STEP) FC2=F(L+1) GC2=G(L+1) D=FC1*GC2-FC2*GC1 AC=(Y1S*GC2-Y2S*GC1)/D BC=(FC1*Y2S-FC2*Y1S)/D PHASE=ATAN2(BC,AC) PHASE=PHASE*180.D0/PI CC=DSQRT(AC*AC+BC*BC) CN=DSQRT(P1/(PI*RK))/CC DO 78 K=1,NMAX R=K*H FUN(K)=FUN(K)*CN/R 78 CONTINUE RETURN END FUNCTION FW(ETA,RM,Z) IMPLICIT REAL*8(A-H,O-Z) DATA EPS/1.D-9/ RM2=RM*RM FW=1.D0 CNT=1.D0 IF(DABS(ETA).LT.1.D-13) GO TO 3 N=1 GO TO 1 2 FW=FW+CNT N=N+1 1 X=ETA+N-0.5D0 X2=X*X CN=(RM2-X2)/(N*Z) IF(DABS(CN).GE.1.D0) RETURN CNT=CNT*CN IF(DABS(CNT).LT.EPS) RETURN GO TO 2 3 K=RM DO 4 N=1,K X=N-0.5D0 X2=X*X CN=(RM2-X2)/(N*Z) CNT=CNT*CN 4 FW=FW+CNT RETURN END c----------------------------------------------------------------- SUBROUTINE RCWFN(RHO,ETA,MINL,MAXL,FC,FCP,GC,GCP,ACCUR,STEP) c------------------------------------------------------------------ IMPLICIT REAL*8(A-H,O-Z) REAL*8 K,K1,K2,K3,K4,M1,M2,M3,M4 DIMENSION FC(1),FCP(1),GC(1),GCP(1) C *** COULOMB WAVEFUNCTIONS C *** CONTINUED-FRACTION METOD OF STEED C *** A.R.BARNETT ET AL. COMP.PHYS.COM.8,377(1974) VANYA=1.D65 PACE=STEP ACC =ACCUR IF(PACE.LT.100.D0) PACE=100.D0 IF(ACC.LT.1.D-15.OR.ACC.GT.1.D-6) ACC=1.D-6 R = RHO KTR = 1 LMAX= MAXL LMIN1= MINL+1 XLL1=DFLOAT(MINL*LMIN1) ETA2= ETA*ETA TURN= ETA+DSQRT(ETA2+XLL1) IF(R.LT.TURN.AND.DABS(ETA).GE.1.D-6) KTR=-1 KTRP = KTR GO TO 2 1 R = TURN TF = F TFP = FP LMAX = MINL KTRP = 1 2 ETAR = ETA*R RHO2 = R*R PL=DFLOAT(LMAX+1) PMX = PL+0.5D0 C *** CONTINUED FRACTION FOR FP(MAXL)/F(MAXL) FP = ETA/PL+PL/R DK = ETAR*2.D0 DEL = 0.D0 D = 0.D0 F = 1.D0 K = (PL*PL-PL+ETAR)*(2.D0*PL-1.D0) IF(PL*PL+PL+ETAR.NE.0.D0) GO TO 3 R=R+1.D-6 GO TO 2 3 H = (PL*PL+ETA2)*(1.D0-PL*PL)*RHO2 K = K+DK+6.D0*PL*PL D = 1.D0/(D*H+K) DEL = DEL*(D*K-1.D0) IF(PL.LT.PMX) DEL=-R*(PL*PL+ETA2)*(PL+1.D0)*D/PL PL = PL+1.D0 FP=FP+DEL IF(D.LT.0.D0) F=-F IF(PL.GT.20000.D0) GO TO 11 IF(DABS(DEL/FP).GE.ACC) GO TO 3 FP=F*FP IF(LMAX.EQ.MINL) GO TO 5 FC(LMAX+1)=F FCP(LMAX+1)=FP C *** DOWNWARD RECURCION TO MINL FOR F AND FP L = LMAX DO 4 LP=LMIN1,LMAX PL=DFLOAT(L) GC(L+1) =ETA/PL+PL/R GCP(L+1) =DSQRT(ETA2+PL*PL)/PL FC(L) =(GC(L+1)*FC(L+1)+FCP(L+1))/GCP(L+1) FCP(L) =GC(L+1)*FC(L)-GCP(L+1)*FC(L+1) 4 L = L-1 F = FC(LMIN1) FP = FCP(LMIN1) 5 IF(KTRP.EQ.-1) GO TO 1 C *** REPEAT FOR R=TURN IF RHO LT TURN C *** C *** NOW OBTAIN P + I*Q FOR MINL FROM CONT. FRACTION (STEED) P = 0.D0 Q = R-ETA PL= 0.D0 AR= -(ETA2+XLL1) AI= ETA BR= 2.D0*Q BI= 2.D0 WI= 2.D0*ETA DR= BR/(BR*BR+BI*BI) DI= -BI/(BR*BR+BI*BI) DP= -(AR*DI+AI*DR) DQ= (AR*DR-AI*DI) 6 P = P+DP Q = Q+DQ PL= PL+2.D0 AR= AR+PL AI= AI+WI BI= BI+2.D0 D = AR*DR-AI*DI+BR DI= AI*DR+AR*DI+BI T = 1.D0/(D*D+DI*DI) DR= T*D DI= -T*DI H = BR*DR-BI*DI-1.D0 K = BI*DR+BR*DI T = DP*H-DQ*K DQ= DP*K+DQ*H DP= T IF(PL.GT.46000.D0) GO TO 11 IF(DABS(DP)+DABS(DQ).GE.(DABS(P)+DABS(Q))*ACC) GO TO 6 P = P/R Q = Q/R C *** SOLVE FOR FP,G,GP,AND NORMALISE F AT L=MINL G = (FP-P*F)/Q GP = P*G-Q*F W = 1.D0/DSQRT(FP*G-F*GP) G = W*G GP = W*GP IF(KTR.EQ.1) GO TO 8 F = TF FP = TFP LMAX = MAXL C *** RUNGE-KUTTA INTEGR. OF G(MINL) AND GP(MINL) C *** INWARD FROM TURN (FOX AND MAYERS 1968 P. 202) IF(RHO.LT.0.2D0*TURN) PACE=999.D0 R3=1.D0/3.D0 H = (RHO-TURN)/(PACE+1.D0) H2 = 0.5D0*H I2=PACE+0.001D0 ETAH=ETA*H H2LL=H2*XLL1 S =(ETAH+H2LL/R)/R-H2 7 RH2 = R+H2 T = (ETAH+H2LL/RH2)/RH2-H2 K1 = H2*GP M1 = S*G K2 = H2*(GP+M1) M2 = T*(G+K1) K3 = H*(GP+M2) M3 = T*(G+K2) M3 = M3+M3 K4 = H2*(GP+M3) RH = R+H S = (ETAH+H2LL/RH)/RH-H2 M4 = S*(G+K3) G = G+(K1+K2+K2+K3+K4)*R3 GP = GP+(M1+M2+M2+M3+M4)*R3 IF(DABS(GP).GT.VANYA.OR.DABS(G).GT.VANYA) GO TO 11 R = RH I2 = I2-1 IF(I2.GE.0) GO TO 7 W = 1.D0/(FP*G-F*GP) C *** UPWARD RECURSION FROM GC(MINL) AND GCP(MINL) C *** RENORMALISE FC,FCP FOR EACH L-VALUE 8 GC(LMIN1)=G GCP(LMIN1)=GP IF(LMAX.EQ.MINL) GO TO 10 DO 9 L=LMIN1,LMAX T = GC(L+1) GC(L+1) = (GC(L)*GC(L+1)-GCP(L))/GCP(L+1) GCP(L+1)= GC(L)*GCP(L+1)-GC(L+1)*T FC(L+1) = W*FC(L+1) 9 FCP(L+1)= W*FCP(L+1) FC(LMIN1) = FC(LMIN1)*W FCP(LMIN1)= FCP(LMIN1)*W RETURN 10 FC(LMIN1) = W*F FCP(LMIN1)= W*FP RETURN 11 W = 0.D0 G = 0.D0 GP= 0.D0 GO TO 8 END SUBROUTINE SINT(VP,N) C RESTRICTIONS: VP(0)=0, N=5 OR >5. IMPLICIT REAL*8(A-H,O-Z) DIMENSION VP(N) DATA KL/1/ N1=N-1 N2=N-2 N3=N-3 IF(KL.EQ.2) GO TO 7 C11=53.D0/360.D0 C10=-11.D0/30.D0 C11M=323.D0/360.D0 C12M=251.D0/720.D0 C12=-19.D0/720.D0 C20=19.D0/30.D0 C21=-37.D0/360.D0 C21M=173.D0/360.D0 C22=11.D0/720.D0 KL=2 7 CONTINUE V1=VP(1) V2=VP(2) V3=VP(3) V4=VP(4) V5=VP(5) VP(1)=C11M*V1+C10*V2+C11*V3+C12*V4 VP(2)=VP(1)+C21M*V1+C20*V2+C21*V3+C22*V4 DO 1 I=3,N2,2 VP(I)=VP(I-1)+C12*V1+C21M*V2+C20*V3+C21*V4+C22*V5 VP(I+1)=VP(I)+C22*V1+C21*V2+C20*V3+C21M*V4+C12*V5 IF(I.EQ.N3) GO TO 2 IF(I.EQ.N2) GO TO 3 V1=V3 V2=V4 V3=V5 V4=VP(I+3) V5=VP(I+4) 1 CONTINUE 2 V1=V2 V2=V3 V3=V4 V4=V5 V5=VP(N) VP(N1)=VP(N2)+C22*V1+C21*V2+C20*V3+C21M*V4+C12*V5 3 VP(N)=VP(N1)+C12*V1+C11*V2+C10*V3+C11M*V4+C12M*V5 RETURN END