IMPLICIT REAL*8(A-H,O-Z) CHARACTER*35 FSPETTRA DIMENSION PKI(80),SPL(2,4,3,80,25),E2F(4,25),NEMA(4), +EMI(4),EMA(4,),SPT(3,80,25),BDE(2),TMK(80,25) COMMON/ENE1/E2F DATA NEMA/24,12,12,12/ DATA EMI/0.D0,12.25D0,50.D0,300.D0/ DATA EMA/12.25D0,50.D0,300.D0,600.D0/ cc DSQRT(AMAS2) = invariant mass of the final state cc AMA= maas of the triton of He3, BDE(NN(2))=binding energy,NN(2)= Z EMAX=DSQRT(AMAS2)-AMA-BDE(NN(2)) EM1=EMAX CALL LEGGI(SPL,EUP,PUP,PKI,NEMA,FSPETTRA) c DO 51 INUC=1,2 DO 28 IT=1,4 NEP1=NEMA(IT)+1 DO 20 IE=1,NEP1 DO 21 IK=1,NKMAX cc IF(E2F(IE).GT.EM1) the extrema of integration on k become equal, and then cc the integration over the spectral function vanishes. cc cc IF(E2F(IE).LE.EM1) one has to evaluate k_max and k_min and then cc TMK(80,25), i.e. the values of nucleon the 3-momentum, where the cc Spectral functions have to be interpolated. TMK(IK,IE)=...... 20 CONTINUE CALL FSKEI(INUC1,NKMAX,IT,NEP1,EM1,SPL,PKI,SPT,TMK,E2F) cc is on the variable k(===> TMK), while the integration over E is the outer cc one. This explains why we choose the Gauss-legendre points for E and cc the Chebichev points (most suitable for interpolation) for PKI. DO 120 IE=1,NEP1 DO 121 IK=1,NKMAX PKE=SPT(3,IK,IE) F2KE=SPT(2,IK,IE) F1KE=SPT(1,IK,IE) 121 CONTINUE 120 CONTINUE 28 CONTINUE 51 CONTINUE SUBROUTINE LEGGI(SPL,EUP,PUP,PKI,NEMA,FSPETTRA) IMPLICIT REAL*8(A-H,O-Z) CHARACTER*35 FSPETTRA DIMENSION PKI(80),SPL(2,4,3,80,25),E2F(4,25),NEMA(4), +EMI(4),EMA(4) COMMON/ENE1/E2F DATA NEMA/24,12,12,12/ DATA EMI/0.D0,12.25D0,50.D0,300.D0/ DATA EMA/12.25D0,50.D0,300.D0,600.D0/ OPEN(9,FILE=FSPETTRA,STATUS='OLD',RECL=110) c A. Kievsky et al Phys.Rev. C56 (1997) 64-75 c EUP= max energy of the spectator pair in MeV c PUP= max momentum of the struck nucleon in 1/fm c NETOT= total number of energies of the spectator pair c NKMAX= total number of momenta of the struck nucleon READ(9,*) EUP,PUP,NETOT,NKMAX c INUC =1 proton spectral functions C INUC =2 neutron spectral functions DO 1 INUC=1,2 NE1=1 IF(INUC.EQ.2) NE1=2 c IT is the number of grids for the enrgies [EMI(IT)-EMA(IT)] DO 1 IT=1,4 NEP1=NEMA(IT)+1 c E2F(IT,IE)=energies in the IT-th grid (suitable for c gaussian integration on that grid) READ(9,*) (E2F(IT,IE),IE=1,NEP1) c PKI(IK)=momenta of the struck nucleon, suitable for the Chebychev c interpolation READ(9,*) (PKI(IK),IK=1,NKMAX) c DO 1 IE=NE1,NEP1 DO 1 IN=1,3 READ(9,*) INN c PKE=SPL(INUC,IT,3,IK,IE)/k unpolarized spectral function c normalized to 4 \pi \int dE/(197.23) \int dk k^2 P(k,E) =1 c F2KE=SPL(INUC,IT,2,IK,IE)/k 0 c F1KE=SPL(INUC,IT,1,IK,IE)/k c PKE,F1KE, F2KE all in fm**4 READ(9,*) E22,(SPL(INUC,IT,IN,IK,IE), IK=1,NKMAX) 1 CONTINUE CLOSE(9) RETURN END SUBROUTINE FSKEI(INUC1,NKMAX,IT,NEP1,EM1,SPL,PKI,SPT,TMK,E2F) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION PKI(80),PK(80),SPL(2,4,3,80,25),SPT(3,80,25),TMK(80,25) +,TP(1),TPN(1),SP1(80,1),SP2(80,1),SP3(80,1),SP1N(80,1),SP2N(80,1) +,SP3N(80,1),E2F(4,25) SAVE DO 1 IE=1,NEP1 DO 1 IK=1,NKMAX DO 1 IN=1,3 SPT(IN,IK,IE)=0. 1 CONTINUE NE1=1 IF(INUC1.EQ.2) NE1=2 DO 3 IE=NE1,NEP1 IF(EM1.LT.E2F(IT,IE)) RETURN DO 5 IK=1,NKMAX PK(IK)=TMK(IK,IE) 5 CONTINUE DO 2 IK=1,80 SP1(IK,1)=SPL(INUC1,IT,1,IK,IE) SP2(IK,1)=SPL(INUC1,IT,2,IK,IE) SP3(IK,1)=SPL(INUC1,IT,3,IK,IE) 2 CONTINUE TP(1)=E2F(IT,IE) TPN(1)=E2F(IT,IE) CALL INTER(1,PKI,TP,PK,TPN,SP1,SP1N,80,1,NKMAX,1) CALL INTER(1,PKI,TP,PK,TPN,SP2,SP2N,80,1,NKMAX,1) CALL INTER(1,PKI,TP,PK,TPN,SP3,SP3N,80,1,NKMAX,1) DO 152 IK=1,NKMAX SPT(1,IK,IE)=SP1N(IK,1)/PK(IK) SPT(2,IK,IE)=SP2N(IK,1)/PK(IK) SPT(3,IK,IE)=SP3N(IK,1)/PK(IK) C WRITE(6,*) SPT(3,IK,IE),IK,IE 152 CONTINUE C IF(INUC1.EQ.2.OR.IT.GT.1) GO TO 3 C DO 190 IN=1,3 C INN=IN C IF(IN.EQ.3) INN=0 C WRITE(8,161) INN C161 FORMAT(2X,' FUNZIONE K*F',I1,//) C WRITE(6,141) (PKI(IK),IK=1,80) C WRITE(6,143) E2F(IT,IE),(SPL(INUC1,IT,IN,IK,IE), IK=1,80) C WRITE(6,141) (PK(IK),IK=1,NKMAX) C WRITE(6,143) E2F(IT,IE),(SPT(IN,IK,IE), IK=1,NKMAX) C 190 CONTINUE 3 CONTINUE C 141 FORMAT(10X,10F12.3,/10X,10F12.3) C143 FORMAT(2X,F8.3,10D12.5,/10X,10D12.5,/10X,10D12.5,/10X,10D12.5,/10X C +,10D12.5) RETURN END SUBROUTINE INTER(KD,XI,YJ,X,Y,FN,F,NN,MMN,NA,MA) C INTERPOLATION OF A FUNCTION DEPENDING ON KD VARIABLES (KD=1,2) C by THE LAGRANGE POLYNOM METHOD (SEE I.S. BEREZIN AND C N.P.ZHIDKOV 'COMPUTING METHODS' , PERGAMON PRESS, P.52 ) C C XI(I), YJ(J) ARE THE COORDINATES WHERE THE FUNCTION IS KNOWN C ( F(XI(I),YJ(J)) = FN(I,J)). IT IS BETTER TO USE THE C CHEBYCHEV POINTS. C NN AND MN ARE THE MAXIMAL NUMEBERS OF XI(IN) AND YJ(JN) C IN=1,2,...,NN ; JN=1,2,...,MN C X(L1) , Y(L2) ARE THE COORDINATES WHERE THE FUNCTION MUST BE INTERPOLATED C (F(X(L1),Y(L2))= F(L1,L2)). C NA AND MA ARE THE MAXIMAL NUMEBERS OF X(L1) AND Y(L2) C L1=1,2,...,NA ; L2=1,2...,MA C C F(L1,L2)=LNM(X(L1),Y(L2))= SUM_{ON I AND J} (FN(I,J)*OME(X)*OME(Y)/ C ((X-XI(I))*(Y-YJ(J))*OME'(XI(I))*OM'(YJ(J))) ) C WITH C OME(X)=(X-XI(1))*(X-XI(2))*...*(X-XI(NN)) C OME(Y)=(Y-YJ(1))*(Y-YJ(2))*...*(Y-YJ(MN)) C AND OME' THE CORRESPONDING DERIVATIVE C IMPLICIT REAL*8(A-H,O-Z) DATA ICALL/0/ DIMENSION X(1),Y(1),XI(1),YJ(1),FN(100,1),F(80,1), 1COEX(100),COEY(100),PY(100) SAVE C ICALL=ICALL+1 N=NA M=MA MN=MMN IF(KD.EQ.1) M=1 IF(KD.EQ.1) MN=1 C IF(ICALL.GT.1) GO TO 200 C CALCOLO DI ON'(XI(K1))=COEX(K1) DO 1 K1=1,NN COEX(K1)=1. DO 12 K1P=1,NN IF(K1.EQ.K1P) GO TO 12 COEX(K1)=(XI(K1)-XI(K1P))*COEX(K1) IF(COEX(K1).EQ.0.) + WRITE(6,999)K1 12 CONTINUE 1 CONTINUE C CALCOLO DI OM'(YJ(K2))=COEY(K2) DO 4 K2=1,MN COEY(K2)=1. IF (KD.EQ.1) GO TO 4 DO 13 K2P=1,MN IF(K2.EQ.K2P) GO TO 13 COEY(K2)=(YJ(K2)-YJ(K2P))*COEY(K2) IF(COEY(K2).EQ.0.) WRITE(6,999)K2 999 FORMAT(2X,'COE= 0 IN INTERP',5X,I10) 13 CONTINUE 2 CONTINUE 4 CONTINUE 200 CONTINUE DO 9 L1=1,N DO 9 L2=1,M FSUM=0. DO 11 K2=1,MN C CALCOLO DI OMJ(Y) PY(K2)=1. IF(KD.EQ.1) GO TO 11 DO 5 K2P=1,MN IF(K2.EQ.K2P) GO TO 5 PY(K2)=PY(K2)*(Y(L2)-YJ(K2P)) 5 CONTINUE 6 CONTINUE PY(K2)=PY(K2)/COEY(K2) 11 CONTINUE DO 7 K1=1,NN C CALCOLO DI ONI(X) PX=1. DO 3 K1P=1,NN IF(K1.EQ.K1P) GO TO 3 PX=PX*(X(L1)-XI(K1P)) 3 CONTINUE PX=PX/COEX(K1) DO 7 K2=1,MN FSUM=FSUM+FN(K1,K2)*PX*PY(K2) 7 CONTINUE F(L1,L2)=FSUM 9 CONTINUE RETURN END