! CEBAF A. Deur 9/01/98 code to generate the input file of positron.f SUBROUTINE POS_WISER(Z,N,E_IN,P_IN,TH_IN,radlen_IN,radlen_after, > XSEC) ! Z, N are target info ! E_IN is beam energy in GeV ! P_IN is scattered particle momentum in GeV ! TH is scattering angle in rad ! radlen_IN is total radiation length in % before reaction, used for estimating pi0 production ! radlen_after is total radiation length in % after reaction, using for estimating pi0 decay ! XSEC is the output cross section in nbarn/GeV/sr, per nucleon ! IMPORTANT: to get the cross section per nuclei, times this by ! (A)**0.8 ! ! X. Zheng Oct 2007 IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER Z,N CHARACTER*3 PART COMMON/QD/QDF COMMON/DEL/IP C DIMENSION C(0:4,8) C DIMENSION A(0:4) real use1(1000),use2(1000) ! Note: mass are used in MeV only in wiser_sub.f and pos_wiser_sub.f DATA PI/3.1415926536/,AMP/938.28/,AMD/1875.63/,AMPI/139.6/ DATA AMPIZ/134.9734/ logical POSDEBUG c POSDEBUG=.true. POSDEBUG=.false. E1=E_IN*1000. P=P_IN*1000. THP=TH_IN*180./PI radlen_1=radlen_IN radlen_2=radlen_after if (POSDEBUG) print *,'in pos_wiser, radlen_before=',radlen_1, > radlen_IN if (POSDEBUG) print *,'in pos_wiser, radlen_after=',radlen_2, > radlen_after epos = P IA=Z+N IF(IA.EQ.2.OR.IA.EQ.3.OR.IA.EQ.4.OR.IA.EQ.12.OR.IA.EQ.16)THEN ELSEIF(IA.EQ.1)THEN PRINT *, ' NUCLEON CASE' ELSE if (POSDEBUG) PRINT 107 107 FORMAT(' NO SPECIFIC SPECTRAL FUNCTION FOR THIS NUCLEUS') if (POSDEBUG) PRINT 108 108 FORMAT(' A = 16 SPECTRAL FUNCTION WILL BE USED') ENDIF C--------------------------------------------------------------------- if (POSDEBUG) print *,'rad_len=',radlen_1 C 'AN' IS EFFECTIVE NUMBER OF NUCLEONS FOR PION PRODUCTION C here we need average of pion+ and pion- to calculate the pion0 C production, regardless of what we are interested before. C--------------------------------------------------------------- c NMAX should be up to pion maximum momentum. Here I take an approximate c Ebeam-500MeV c DELP is the bin size in pion crss section arrays. Here defined as 100MeV DELP=100. NMAX=INT((E1-500.-P)/DELP+1) if (POSDEBUG) print *,'NMAX=',NMAX P=P-DELP do I=1,NMAX P=P+DELP PART='PI+' IF(PART.EQ.'P')THEN AN=N/3.+2.*Z/3. IP=1 ITYPE=5 ELSEIF(PART.EQ.'N')THEN AN=Z/3.+2.*N/3. IP=-1 ELSEIF(PART.EQ.'PI+')THEN AN=Z/3. IP=2 ITYPE=1 ELSEIF(PART.EQ.'PI-')THEN AN=N/3. IP=2 ITYPE=2 ELSEIF(PART.EQ.'PI0')THEN AN=2.*(N+Z)/3. IP=2 ELSE STOP ENDIF if (POSDEBUG) print *,'POSWISER:',PART,AN,IP,ITYPE IF(IA.EQ.1)THEN if (POSDEBUG) PRINT 114 114 FORMAT(' ONLY DELTA MECHANISM INCLUDED') ELSEIF(ABS(IP).EQ.1)THEN IF(IA.GT.1.AND.IA.LT.5)THEN DLF=IA ELSE DLF=7. ENDIF if (POSDEBUG) PRINT 102 102 FORMAT(' ENTER LEVINGER FACTOR OR 0 FOR DEFAULT VALUE') READ *,AL IF(AL.EQ.0.)AL=DLF QDF=AL*N*Z/FLOAT(IA) ELSE ENDIF PTP=P IF(ABS(IP).EQ.1)THEN E=SQRT(PTP**2+AMP**2) TP=PTP**2/(E+AMP) AJ=PTP/E D2QD=0. D2QF=0. ELSE E=SQRT(PTP**2+AMPI**2) TP=PTP**2/(E+AMPI) AJ=PTP/E D2QD=0. D2QF=0. ENDIF D2DEL=0. IF(ABS(IP).EQ.2.OR.IP.EQ.0)THEN Call WISER_ALL_SIG(E1,PTP,THP,RADLEN_1,ITYPE,TOTAL) ! TOTAL is pion cross section in nbarn/GeV/sr for a proton target TOTAL=TOTAL*AN ENDIF if (POSDEBUG) print *,'PI+:',TOTAL use1(I)=PTP use2(I)=TOTAL icount=icount+1 PART='PI-' IF(PART.EQ.'P')THEN AN=N/3.+2.*Z/3. IP=1 ITYPE=5 ELSEIF(PART.EQ.'N')THEN AN=Z/3.+2.*N/3. IP=-1 ELSEIF(PART.EQ.'PI+')THEN AN=Z/3. IP=2 ITYPE=1 ELSEIF(PART.EQ.'PI-')THEN AN=N/3. IP=2 ITYPE=2 ELSEIF(PART.EQ.'PI0')THEN AN=2.*(N+Z)/3. IP=2 ELSE STOP ENDIF if (POSDEBUG) print *,'POSWISER:',PART,AN,IP,ITYPE IF(IA.EQ.1)THEN PRINT 114 ELSEIF(ABS(IP).EQ.1)THEN IF(IA.GT.1.AND.IA.LT.5)THEN DLF=IA ELSE DLF=7. ENDIF PRINT 102 READ *,AL IF(AL.EQ.0.)AL=DLF QDF=AL*N*Z/FLOAT(IA) ELSE ENDIF PTP=P IF(ABS(IP).EQ.1)THEN E=SQRT(PTP**2+AMP**2) TP=PTP**2/(E+AMP) AJ=PTP/E D2QD=0. D2QF=0. ELSE E=SQRT(PTP**2+AMPI**2) TP=PTP**2/(E+AMPI) AJ=PTP/E D2QD=0. D2QF=0. ENDIF D2DEL=0. IF(ABS(IP).EQ.2.OR.IP.EQ.0)THEN Call WISER_ALL_SIG(E1,PTP,THP,RADLEN_1,ITYPE,TOTAL) ! TOTAL is pion cross section in nbarn/GeV/sr for a proton target TOTAL=TOTAL*AN ENDIF if (POSDEBUG) print *,'PI-:',TOTAL C--------------------------------------------------------------- if (POSDEBUG) print *,'use1=',use1(I),' use2=',use2(I),I use1(I)=PTP use2(I)=use2(I)+TOTAL use2(I)=use2(I)/2 if (use2(I).lt.0) use2(I)=0. if (POSDEBUG) print *,'PI0:',use2(I) ENDDO ! END 1,NMAX !radiation length before interaction ( in %) : eff=radlen_1 EPHOTON=AMPIZ/2. IF (POSDEBUG) print *,'POS_WISER_SUB: eff=',eff,' z=',z,' ia=',ia C THICK=0.01995 ! = 5*10^21 atom/cm^2 C EFF=EFF*(THICK/.8304+0.12) ! ONLY ONE BIN FOR EPRIME I=1 RAT=EPOS/EPHOTON CROSS=0. if (POSDEBUG) print *,'NMAX=',NMAX do 90 II=1,NMAX P=use1(II) CRS=use2(II) if (POSDEBUG) print *, II,p,epos if(p.lt.epos) then if (POSDEBUG) print *,'pi0 momentum=',p,' epos=',epos goto 90 ! do not contribute if pi0 momentum is < Eprime endif if (CRS.lt.0) goto 90 EPI=SQRT(AMPIZ**2+P**2) GARMA=EPI/AMPIZ BETA=P/EPI THETA=RTX(GARMA,BETA,rat) C=AING(GARMA,BETA,THETA) c CROSS=CROSS+CRS*C CROSS=CROSS+CRS*C*DELP/1000. ! DELP is pion0 bin size (in MeV) ! used to creat array use1 and use2 ! CRS is in nbarn/GeV/sr, so should divide by 1000 if (POSDEBUG) print *, p,theta,c,crs,cross 90 continue c CROSS=CROSS*0.0446/PI/EPHOTON*1000. c I think 0.0446 is the radiation length after scattering. For 20cm 3He cell c used in A1n (Gore) it is about 0.05 (absolute, not %) if (POSDEBUG) print *,'CROSS=',CROSS, > ' radlen_after=',radlen_after,PI,EPHOTON CROSS=CROSS*radlen_after/PI/EPHOTON*1000. cross=eff*cross ! cross in nanobarn/GeV*str already XSEC=cross 999 continue RETURN END *VTP SUBROUTINE VTP(AMT,AM1,EI,W0,TP,TH,GN) C TIATOR-WRIGHT VIRTUAL PHOTON SPECTRUM C PHYS. REV. C26,2349(1982) AND NUC. PHYS. A379,407(1982) IMPLICIT REAL*8 (A-H,O-Z) DATA AME/.511/,PI/3.1416/ EF0=EI-W0 AKI=SQRT(EI**2-AME**2) AKF0=SQRT(EF0**2-AME**2) AKP=SQRT(TP**2+2.*AM1*TP) EP=TP+AM1 AR=EI+AMT-EP BR=EF0*(AKP*COS(TH)-AKI)/AKF0 BRP=(AKF0/EF0)**2*BR A=AME**2-EI*EF0 B=AKI*AKF0 D=-AME**2*BR*(EI/EF0-1.)/AR AP=A-D BP=B+D AN1=1./137./2./PI*W0**2/AKI**2 APB=-AME**2*(AKI-AKF0)**2/(AME**2+EI*EF0+AKI*AKF0) AN1=AN1*B/BP*(AR+BR)/(AR-AP/BP*BR) AN2=1.-2.*A/W0**2 AN4=((AP-BP)*(AR+BR)/APB/(AR-BR)) IF(AN4.LE.0.)GO TO 1 AN2=AN2*LOG(AN4) AN3=-4.*B/W0**2 ANE=AN1*(AN2+AN3) D0=AMT+EI-EP+EF0/AKF0*(AKP*COS(TH)-AKI) R=(AMT+W0-EP/AKP*W0*COS(TH))/D0 GN=ANE*R/W0 IF(GN.LT.0.)GN=0. RETURN 1 GN=0. RETURN END *DEP SUBROUTINE DEP(E1,TP,TH,IP,D2N) C QUASI-DEUTERON CROSS SECTION IMPLICIT REAL*8 (A-H,O-Z) COMMON/QD/QDF C DIMENSION C(0:4,8),A(0:4) DATA PI/3.1416/,AM/939./,AMD/1876./ PN=SQRT(TP**2+2.*AM*TP) EP=TP+AM CALL KINE(AMD,AM,AM,PN,TH,W0,THC) IF(W0.GE.E1)GO TO 1 IF(W0.LE.0.)GO TO 1 W0G=W0/1000. CALL SIGD(W0G,THC,IP,DSQD) CALL PART(AMD,AM,AM,PN,TH,AJT,AJW) DSQD=AJT*DSQD C CROSS SECTION IN UB/MEV-SR CALL VTP(AMD,AM,E1,W0,TP,TH,PHI) D2N=QDF*PHI*DSQD RETURN 1 D2N=0. RETURN END *SIGD SUBROUTINE SIGD(E,TH,IP,DSQD) C DEUTERON CROSS SECTION C BASED ON FIT OF THORLACIUS & FEARING C PHYS. REV. C33,1830(1986) C ENERGY RANGE 10 - 625 MEV C C E[GEV] IN LAB SYSTEM C TH[RAD] & DSQD[UB/SR] IN CENTER-OF-MOMENTUM SYSTEM C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION C0(8),C1(4),C2(4),C3(4),C4(4) DIMENSION A(0:4),B(4,4) DATA C0/2.61E2,-1.10E2,2.46E1,-1.71E1,5.76E0,-2.05E0,2.67E-1, 1 1.13E2/ DATA C1/1.68E1,-4.66E1,2.56E0,-4.72E0/ DATA C2/-2.03E2,-8.12E1,-4.05E0,-5.99E0/ DATA C3/-1.77E1,-3.74E1,-5.07E-1,-5.40E0/ DATA C4/-2.05E0,-7.05E0,9.40E-1,-2.05E0/ X=COS(TH) IF(E.LE.625.)THEN C TEST FOR NEUTRON X=IP*X C COEFICIENTS A(0)=C0(1)*EXP(C0(2)*E)+ C0(3)*EXP(C0(4)*E) A(0)=A(0)+(C0(5)+C0(6)*E)/(1.+C0(8)*(E-C0(7))**2) DSQD=A(0)*P(0,X) DO 2 L=1,4 B(1,L)=C1(L) B(2,L)=C2(L) B(3,L)=C3(L) 2 B(4,L)=C4(L) DO 1 L=1,4 A(L)=B(L,1)*EXP(B(L,2)*E)+ B(L,3)*EXP(B(L,4)*E) 1 DSQD=DSQD+A(L)*P(L,X) ELSEIF(E.LT.700.)THEN DSDQ=.3 ELSEIF(E.LT.800.)THEN DSDQ=.15 ELSEIF(E.LT.900.)THEN DSDQ=.1 ELSE DSDQ=55./(E-350.) ENDIF RETURN END *LEG REAL*8 FUNCTION P(L,X) C LEGENDRE POLYNOMIALS IMPLICIT REAL*8 (A-H,O-Z) IF(L.EQ.0)THEN P=1. ELSEIF(L.EQ.1)THEN P=X ELSEIF(L.EQ.2)THEN P=.5*(3.*X**2-1.) ELSEIF(L.EQ.3)THEN P=.5*(5.*X**3-3.*X) ELSEIF(L.EQ.4)THEN P=1./8.*(35.*X**4-30.*X**2+3.) ELSE P=0. ENDIF RETURN END *DELTA SUBROUTINE DELTA(E1,TP,TH,D2DEL) C PHOTOPRODUCTION OF NUCLEONS AND PIONS VIA DELTA IMPLICIT REAL*8 (A-H,O-Z) COMMON/DEL/IP DATA PI/3.1416/,AM/939./,AMPI/139./ IF(ABS(IP).EQ.1)THEN AM1=AM AM2=AMPI ELSE AM1=AMPI AM2=AM ENDIF EP=TP+AM1 PN=SQRT(EP**2-AM1**2) CALL KINE(AM,AM1,AM2,PN,TH,W,TC) IF(W.LE.0.)GO TO 1 IF(W.GE.E1)GO TO 1 CALL PART(AM,AM1,AM2,PN,TH,AJT,AJW) CALL SIGMA(W,TC,DSIGG) CALL VTP(AM,AM1,E1,W,TP,TH,PHI) D2DEL=PHI*DSIGG*AJT C CROSS SECTION IN UB/MEV-SR RETURN 1 D2DEL=0. RETURN END *PART SUBROUTINE PART(AMT,AM1,AM2,PN,TN,AJT,AJW) IMPLICIT REAL*8 (A-H,O-Z) C PARTIAL DERIVATIVES DATA PI/3.1416/ DT=PI/50. DP=10. C ANGLE TNP=TN+DT TNM=TN-DT CALL KINE(AMT,AM1,AM2,PN,TNP,W,TCP) CALL KINE(AMT,AM1,AM2,PN,TNM,W,TCM) AJT=(COS(TCP)-COS(TCM))/(COS(TNP)-COS(TNM)) AJT=ABS(AJT) C ENERGY PNP=PN+DP PNM=PN-DP CALL KINE(AMT,AM1,AM2,PNP,TN,WP,TC) CALL KINE(AMT,AM1,AM2,PNM,TN,WM,TC) AJW=(WP-WM)/(PNP-PNM) AJW=ABS(AJW) RETURN END *KINE SUBROUTINE KINE(AMT,AM1,AM2,PN,TH,W,TC) C COMPUTES CM VARIABLES FROM LAB VARIABLES IMPLICIT REAL*8 (A-H,O-Z) EP=SQRT(PN**2+AM1**2) PNT=PN*SIN(TH) PNL=PN*COS(TH) ANUM=PN**2+AM2**2-(AMT-EP)**2 DEN=2.*(PNL+AMT-EP) W=ANUM/DEN IF(W.LE.0.)W=0. C INVARIANT MASS WW=SQRT(AMT**2+2.*W*AMT) C CM VARIABLES PCT=PNT B=W/(AMT+W) G=(W+AMT)/WW PCL=G*(PNL-B*EP) PCS=PCL**2+PCT**2 PC=SQRT(PCS) CTHC=PCL/PC TC=ACOS(CTHC) RETURN END *SIGMA SUBROUTINE SIGMA(E,THRCM,SIGCM) IMPLICIT REAL*8 (A-H,O-Z) GAM=100. PI=ACOS(-1.) IF(E.GT.420.)THEN SIGCM=(1.+420./E)*90./4./PI ELSE SIGCM=360.*(5.-3.*COS(THRCM)**2)/16./PI/(1.+(E-320.)**2/GAM**2) ENDIF RETURN END C BEGIN QUASI-FREE SCATTERING CROSS SECTION *EP SUBROUTINE EP(E1,TP,THP,DSEP) C ELECTRO PROTON PRODUCTION CROSS SECTIONS IMPLICIT REAL*8 (A-H,O-Z) COMMON/Pcm/PH(10),WPH(10) DATA AML/.511/,PI/3.14159265/ CALL GAUSAB(10,PH,WPH,0.D0,2.*PI,PI) AK=SQRT(E1**2-AML**2) CALL SEP(AK,TP,THP,DSEP) DSEP=DSEP*1.E4 C CROSS SECTION IN UB/MEV-SR END *DOT REAL*8 FUNCTION DOT(V,U) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION V(3),U(3) DOT=0. DO 1 I=1,3 1 DOT=DOT+V(I)*U(I) RETURN END *CROSS SUBROUTINE CROSS(V,U,W) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION V(3),U(3),W(3) W(1)=V(2)*U(3)-V(3)*U(2) W(2)=V(3)*U(1)-V(1)*U(3) W(3)=V(1)*U(2)-V(2)*U(1) RETURN END *GAUSAB C SUBROUTINE GAUSAB C SUBROUTINE GAUSAB(N,E,W,A,B,C) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION E(24),W(24) DATA PI/3.141592653589793238462643D0/,EPS/1.D-16/ IF(A.GE.C.OR.C.GE.B)STOP C STOPS PROGRAM IF A, C, B ARE OUT OF SEQUENCE AL=(C*(A+B)-2*A*B)/(B-A) BE=(A+B-2*C)/(B-A) M=(N+1)/2 DN=N DO 5 I=1,M DI=I X=PI*(4.D0*(DN-DI)+3.D0)/(4.D0*DN+2.D0) XN=(1.D0-(DN-1.D0)/(8.D0*DN*DN*DN))*COS(X) IF(I.GT.N/2) XN=0 DO 3 ITER=1,10 X=XN Y1=1.D0 Y=X IF(N.LT.2) GO TO 2 DO 1 J=2,N DJ=J Y2=Y1 Y1=Y 1 Y=((2.D0*DJ-1.D0)*X*Y1-(DJ-1.D0)*Y2)/DJ 2 CONTINUE YS=DN*(X*Y-Y1)/(X*X-1.D0) H=-Y/YS XN=X+H IF(ABS(H).LT.EPS) GO TO 4 3 CONTINUE 4 E(I)=(C+AL*X)/(1.D0-BE*X) E(N-I+1)=(C-AL*X)/(1.D0+BE*X) GEW=2.D0/((1.D0-X*X)*YS*YS) W(I)=GEW*(AL+BE*C)/(1.D0-BE*X)**2 W(N-I+1)=GEW*(AL+BE*C)/(1.D0+BE*X)**2 5 CONTINUE RETURN END *VECT SUBROUTINE VECT(THP,THE,PHI,P,AK1,AK2) C CARTESIAN COMPONENTS OF ELECTRON AND PROTON VECTORS IMPLICIT REAL*8 (A-H,O-Z) COMMON/V/AK1V(3),AK2V(3),QV(3),PV(3),PP(3) PV(1)=P*SIN(THP) PV(2)=0. PV(3)=P*COS(THP) AK1V(1)=0. AK1V(2)=0. AK1V(3)=AK1 AK2V(1)=AK2*SIN(THE)*COS(PHI) AK2V(2)=AK2*SIN(THE)*SIN(PHI) AK2V(3)=AK2*COS(THE) QV(1)=AK1V(1)-AK2V(1) QV(2)=AK1V(2)-AK2V(2) QV(3)=AK1V(3)-AK2V(3) PP(1)=PV(1)-QV(1) PP(2)=PV(2)-QV(2) PP(3)=PV(3)-QV(3) RETURN END *AMAG REAL*8 FUNCTION AMAG(V) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION V(3) AMAG=0. DO 1 I=1,3 1 AMAG=AMAG+V(I)**2 AMAG=SQRT(AMAG) RETURN END *LEPT SUBROUTINE LEPT(E1,E2,AK1,AK2,AML,QS,QUS,THE,V) C LEPTON FACTORS FOR COINCIDENCE CROSS SECTION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION V(5) V(1)=(QUS/QS)**2*(E1*E2+AK1*AK2*COS(THE)+AML**2) X=AK1*AK2*SIN(THE) V(2)=X**2/QS+QUS/2. V(3)=QUS/QS*X/SQRT(QS)*(E1+E2) V(4)=X**2/QS V(5)=0. RETURN END *D4S SUBROUTINE D4S(AK1,AK2,THE,P,PP,THQP,CPHIP,DSIG) C FULLY DIFFERENTIAL CROSS SECTION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION V(5),W(5) DATA AM/939./,AML/.511/,PI/3.14159265/,A/855./ QS=AK1**2+AK2**2-2.*AK1*AK2*COS(THE) E1=SQRT(AK1**2+AML**2) E2=SQRT(AK2**2+AML**2) QUS=2.*(E1*E2-AK1*AK2*COS(THE)-AML**2) SM=2.*(1.44)**2/QUS**2*AK2/AK1 EP=SQRT(AM**2+P**2) PS=EP*P FNS=1./(1.+QUS/A**2)**4 CALL LEPT(E1,E2,AK1,AK2,AML,QS,QUS,THE,V) CALL FORM(QS,QUS,P,THQP,CPHIP,W) SUM=0. DO 1 I=1,5 1 SUM=SUM+V(I)*W(I) DSIG=SM*PS*FNS*SUM*SGSL(PP) RETURN END *STHE SUBROUTINE STHE(D2S) C INTEGRAL OVER ELECTRON POLAR ANGLE IMPLICIT REAL*8 (A-H,O-Z) COMMON/S/ AK1,AK2,THE,P,THP COMMON/E/TH1(12),WT1(12),TH2(12),WT2(12) COMMON/E1/TH3(24),WT3(24) D2S1=0. DO 1 I=1,12 THE=TH1(I) CALL SPHI(D3S) 1 D2S1=D2S1+D3S*WT1(I)*SIN(THE) D2S2=0. DO 2 I=1,12 THE=TH2(I) CALL SPHI(D3S) 2 D2S2=D2S2+D3S*WT2(I)*SIN(THE) D2S3=0. DO 3 I=1,24 THE=TH3(I) CALL SPHI(D3S) 3 D2S3=D2S3+D3S*WT3(I)*SIN(THE) D2S=D2S1+D2S2+D2S3 RETURN END *SPHI SUBROUTINE SPHI(D3S) C INTEGRATE OVER ELECTRON AZIMUTHAL ANGLE IMPLICIT REAL*8 (A-H,O-Z) COMMON/S/ AK1,AK2,THE,P,THP COMMON/V/AK1V(3),AK2V(3),QV(3),PV(3),PP(3) COMMON/Pcm/PH(10),WPH(10) DIMENSION QXP(3),AK1X2(3) D3S=0. DO 1 I=1,10 PHI=PH(I) CALL VECT(THP,THE,PHI,P,AK1,AK2) CALL CROSS(QV,PV,QXP) CALL CROSS(AK1V,AK2V,AK1X2) C PROTON THETA CTHEP=DOT(PV,QV)/AMAG(PV)/AMAG(QV) THQP=ACOS(CTHEP) C PROTON PHI CPHIP=DOT(QXP,AK1X2)/AMAG(QXP)/AMAG(AK1X2) PPM=AMAG(PP) CALL D4S(AK1,AK2,THE,P,PPM,THQP,CPHIP,DSIG) 1 D3S=D3S+DSIG*WPH(I) RETURN END *FORM SUBROUTINE FORM(QS,QUS,P,THQP,CPHIP,W) C NUCLEAR FORM FACTORS IMPLICIT REAL*8 (A-H,O-Z) COMMON/M/ZZ,NN COMMON/DEL/IP DIMENSION W(5) DATA AM/939./,UP/2.79/,UN/-1.91/ IF(IP.EQ.1)THEN Z=ZZ N=0. ELSEIF(IP.EQ.-1)THEN Z=0. N=NN ELSE Z=0. N=0. ENDIF Y=P/AM*SIN(THQP) W(1)=Z W(2)=Z*Y**2 W(2)=W(2)+(Z*UP**2+N*UN**2)*QS/2./AM**2 W(3)=-2.*Z*Y*CPHIP W(4)=Z*Y**2*(2.*CPHIP**2-1.) W(5)=0. RETURN END *SEP SUBROUTINE SEP(AK,TP,THPP,D2S) IMPLICIT REAL*8 (A-H,O-Z) COMMON/S/AK1,AK2,THE,P,THP COMMON/E/TH1(12),WT1(12),TH2(12),WT2(12) COMMON/E1/TH3(24),WT3(24) DATA AM/939./,AML/.511/,BE/16./,PI/3.1416/ THP=THPP AK1=AK AK2=AK1-TP-BE C GAUSSIAN POINTS FOR THE THEMAX=AML*(AK1-AK2)/AK1/AK2 CALL GAUSAB(12,TH1,WT1,0.D0,2.*THEMAX,THEMAX) CALL GAUSAB(12,TH2,WT2,2.*THEMAX,100.*THEMAX,10.*THEMAX) A3=100.*THEMAX C3=A3+(PI-A3)/10. CALL GAUSAB(24,TH3,WT3,A3,PI,A3+(PI-A3)/10.) P=SQRT(TP**2+2.*AM*TP) CALL STHE(D2S) IF(AK2.LE.0.)D2S=0. RETURN END *SGSL REAL*8 FUNCTION SGSL(P) C P INTEGRAL OVER SGSL NORMALIZED TO 1/4PI IMPLICIT REAL*8 (A-H,O-Z) COMMON/SP/IA IF(IA.EQ.2)THEN C BEGIN 2-H PP=P/197.3 SGS=3.697-7.428*PP-2.257*PP**2 SGS=SGS+3.618*PP**3-1.377*PP**4+.221*PP**5-.013*PP**6 IF(SGS.LT.-293.)GO TO 1 SGS=EXP(SGS) SGS=SGS/.18825/4./3.1416/(197.3)**3 SGSL=SGS/1. ELSEIF(IA.EQ.3)THEN C BEGIN 3-HE IF(-(P/33)**2.LT.-293.)GO TO 1 SGS=2.4101E-6*EXP(-P/33) SGS=SGS-1.4461E-6*EXP(-(P/33)**2) SGS=SGS+1.6871E-10*EXP(-(P/493)**2) SGSL=SGS/2. ELSEIF(IA.EQ.4)THEN C BEGIN 4-HE IF(-(P/113.24)**2.LT.-293.)GO TO 1 SGS=1.39066E-6*EXP(-(P/113.24)**2) SGS=SGS+3.96476E-9*EXP(-(P/390.75)**2) SGSL=SGS/2. SGSL=SGSL/2./3.1416 ELSEIF(IA.EQ.12)THEN C BEGIN 12-C IF(-(P/127)**2.LT.-293.)GO TO 1 SGS=1.7052E-7*(1.+(P/127)**2)*EXP(-(P/127)**2) SGS=SGS+1.7052E-9*EXP(-(P/493)**2) SGSL=SGS/6. ELSE C BEGIN 16-O IF(-(P/120)**2.LT.-293.)GO TO 1 SGS=3.0124E-7*(1.+(P/120)**2)*EXP(-(P/120)**2) SGS=SGS+1.1296E-9*EXP(-(P/493)**2) SGSL=SGS/8. ENDIF RETURN 1 SGSL=0. RETURN END C ----------- from JPChen's positron code------------- FUNCTION AING(GM,BT,TH) IMPLICIT REAL*8 (A-H,O-Z) DT=TH/500. AING=0. THE=-DT/2. DO I=1,500 THE=THE+DT Y=1./FUNC(GM,BT,THE) AING=AING+Y*DT ENDDO RETURN END FUNCTION FUNC(A,B,THE) IMPLICIT REAL*8 (A-H,O-Z) FUNC=SQRT(SIN(THE)**2+A**2*(COS(THE)+B)**2) RETURN END FUNCTION RTX(A,B,rat) IMPLICIT REAL*8 (A-H,O-Z) PI=3.1415926536 RTX=0. XA=0. XB=PI FA=RAT-FUNC(A,B,XA) IF(ABS(FA).LT.1.E-3)THEN RTX=XA RETURN ENDIF FB=RAT-FUNC(A,B,XB) IF(ABS(FB).LT.1.E-3)THEN RTX=XB RETURN ENDIF IF(FA*FB.GE.0.0)RETURN 100 X=XA-(XB-XA)/(FB-FA)*FA FX=RAT-FUNC(A,B,X) IF(ABS(FX).LT.1.E-3)THEN RTX=X RETURN ENDIF IF(FA*FX.LT.0.)THEN XB=X FB=FX GOTO 100 ELSE XA=X FA=FX GOTO 100 ENDIF END C note: need subroutines from wiser_sub.f