c c c include some modifications for the (e,e') case... c csinpmc234567 program pwba implicit real*8 (a-h,o-z) c real*8 step,r,hst,fac,s2(-1:1),fff(3),fwc(3) character*20 rpiv character*60 ipivot character*26 infile,outfile,fil1,fil2,fil3,fil4,fil5,fil6 PI=dacos(-1.d0) ALPHA=1.d0/137.036d0 hc=197.33d0 c print*,'mete escale' c read(*,*)esc c hc=hc/esc amassp=938.3d0/hc c amassp=0.1d0 alph=137.d0 ua=931.4944/hc sqr2=dsqrt(2.d0) sqr6=dsqrt(6.d0) fourpi=4.d0*pi c call cputim(it1) cpmfac=1.d0 sin2wb=0.227 c sin2wb=0.2330 CV=4*sin2wb-1.d0 CA=1.d0 open(1,file='pvfermi2.in',status='old') 1234 format(60A) 123 format(26A) c call cputim(itbeg) c c c Entrada de datos c c 9991 format(7e13.5) read(1,1234)ipivot read(1,1234)ipivot read(1,123)outfile open(13,file='rfg.out') read(1,1234)ipivot write(13,*)ipivot read(1,1234)ipivot write(13,*)ipivot read(1,*)z1num,z1mass,z2num,z2mass write(13,9991)z1num,z1mass,z2num,z2mass read(1,1234)ipivot write(13,*)ipivot read(1,*)nb,lb2,jb2,spft,iff,ireac,icharg,asig if(ireac.eq.0)then if(icharg.eq.1)then rpiv=' (e,e''p) ' else rpiv=' (e,e''n) ' endif else if(icharg.eq.1) then if(asig.gt.0.)then rpiv=' (nu,nu''p) ' else rpiv=' (anu,anu''p) ' endif else if(asig.gt.0.)then rpiv=' (nu,nu''n) ' else rpiv=' (anu,anu''n) ' endif endif endif c c c write(13,*)nb,lb2,jb2,spft,iff,ireac,icharg,asig read(1,1234)ipivot write(13,*)ipivot read(1,*)Tpfin,tpffin,tpfstep,qinp,pmfin,dpm write(13,9991)tpfin,tpffin,tpfstep,qinp,pmfin,dpm wmini=wmini/hc read(1,1234)ipivot write(13,*)ipivot read(1,*)eki,Emiss,tick write(13,9991)eki,Emiss,tick read(1,1234)ipivot write(13,*)ipivot read(1,123)infile write(13,123)infile read(1,1234)ipivot write(13,*)ipivot read(1,1234)ipivot write(13,*)ipivot read(1,*)step,rfin,nstcos,cosmin,cosmax,imp write(13,1111)step,rfin,nstcos,cosmin,cosmax,imp 1111 format(2e13.6,i5.3,2e13.5,i5.3) read(1,1234)ipivot write(13,*)ipivot read(1,123)fil1 read(1,123)fil2 c read(1,123)fil3 write(13,*)fil1 write(13,*)fil2 c write(13,*)fil3 read(1,1234)ipivot write(13,*)ipivot read(1,123)fil4 read(1,123)fil5 read(1,123)fil6 write(13,*)fil4 write(13,*)fil5 write(13,*)fil6 read(1,1234)ipivot write(13,*)ipivot c read(1,1234)ipivot c write(13,*)ipivot read(1,123)fil3 write(13,*)fil3 read(1,1234)ipivot write(13,*)ipivot read(1,1234)ipivot write(13,*)ipivot read(1,*)ipl,ipt,icc,iwoper,f2s,g1s write(13,*)ipl,ipt,icc,iwoper ,f2s,g1s c c c Fin de la entrada de datos c open(26,file='sigcc1.out') c open(11,file=infile,status='old') open(37,file='proton.out') open(38,file='neutron.out') open(36,file='0int.out') open(39,file='inclus.out') open(57,file='formf.out') open(40,file='sigma.out') open(50,file='pcresp.out') open(51,file='pvresp.out') write(*,*)' ---------------------------------------' write(*,*)' - Program RFG- INCLUSIV- J.M. UDIAS -' write(*,*)' -AMSTERDAM -1994---- Current v. 2.01-' write(*,*)' ---------------------------------------' write(39,*)' ---------------------------------------' write(39,*)' - Program RFG- INCLUSIV- J.M. UDIAS -' write(39,*)' -AMSTERDAM -1994---- Current v. 2.01-' write(39,*)' ---------------------------------------' write(39,*)'-----------------------------------------' write(39,*)' Fermi Momentum=',pmfin,' MeV ' write(39,*)' Binding Energy=',b ,' MeV ' write(39,*)' energy incident=',eki,' MeV' write(39,*)' values of response functions given for ' write(39,*)' w=',tpfin,' MeV' write(39,*)' to w=',tpffin, ' MeV ' write(39,*)'(or the values constrained by kinematics)' write(39,*)' in steps of ',tpfstep,' MeV ' write(39,*)' Ef integration done in ',nstcos, ' steps ' write(39,*)'------------------------------------------' write(*,*)'-----------------------------------------' write(*,*)' Fermi Momentum=',pmfin,' MeV ' write(*,*)' Binding Energy=',b ,' MeV ' write(*,*)' energy incident=', eki, ' MeV ' write(*,*)' values of response functions given for ' write(*,*)' w=',tpfin,' MeV' write(*,*)' to w=',tpffin, ' MeV ' write(*,*)'(or the values constrained by kinematics)' write(*,*)' in steps of ',tpfstep,' MeV ' write(*,*)' Ef integration done in ',nstcos, ' steps ' write(*,*)'------------------------------------------' write(57,*) '--------------------------------------' write(57,*)' Factores de Forma para el caso',rpiv idip=1 if(iff.lt.0)then write(57,*)'dipole dependence set to 1' idip=0 iff=iabs(iff) endif if(iff.eq.2)write(57,*)' F1 set to 0 in the calculation' if(iff.eq.3)write(57,*)' F2 set to 0 in the calculation' if(iff.eq.4)write(57,*)' F1+kF2 set to 0 in the calc.' if(iff.eq.5)write(57,*)' F1=F2=0 only G1=1 enters' write(57,*)'---------------------------------------' write(57,*)' qnu2 (fm-2) ff1 ff2 G1 fcoup ', , ' gm/ge test ' write(57,*)'-----------------------------------------------------' write(38,*) '--------------------------------------' write(38,*)' Funciones de Respuesta para el caso ',rpiv write(38,*)'-----------------------------------------------------' write(38,*)' Pm (MeV) ispf2 mb2 RL RT RLT', , ' RTT RATT ' write(38,*)'-----------------------------------------------------' write(38,*) '--------------------------------------' 4799 format(' Pm spf2 mb2 comp gam sig ', ,' five ') 4759 format('-----------------------------------------------------') write(26,*)' cross section sig_cc1 of de Forest in fm2 ', , rpiv write(26,*)' pm gamma phi theta cc1', , ' L T LT' write(26,*)'-------------------------------------------------' if(iff.eq.4)print*,'F_1 + k F_2 = 0' if(iff.eq.5)print*,'-k F_2 /2/M = 0' c c masa nucleo residual y blanco c c c energia total, momento y energia cinetica del proton final c emiss=emiss/hc B=emiss q=eki/hc pmfmin=pmfin/hc c c fermi momentum in fm-1 c pfermi=pmfmin c print*,'pfermi=',pfermi*hc c c normalization factor for fermi gas c ffg0=3.d0/4.d0/pi/pfermi**3 q2=q*q do 12345 wloop=tpfin,tpffin,tpfstep w=wloop/hc w2=w*w ang=dabs(qinp) cosnu=dcos(ang*pi/180.d0) ene1=(2*w*(1.d0-cosnu)+dsqrt(4.*w**2.*(1.-cosnu)**2.- 1 8.*(1-cosnu)*(w2-q2)))/4./(1-cosnu) ene2=(2*w*(1.d0-cosnu)-dsqrt(4.*w**2.*(1.-cosnu)**2.- 1 8.*(1-cosnu)*(w2-q2)))/4./(1-cosnu) print*,'ene1,ene2',ene1,ene2 if(ene1.gt.0.)then eki=ene1 else eki=ene2 endif ekf=eki-w ekf2=ekf*ekf c print*,w aa=2*(b-w) bb=(b-w)**2-q*q root=4*aa*aa*bb*bb+4*(4*q*q-aa*aa)*(bb*bb+4*q*q*amassp**2) root=dsqrt(root) den=8*q*q-2*aa*aa e1=2*aa*bb+root e2=2*aa*bb-root e1=e1/den e2=e2/den c root=4*aa*aa*bb*bb+4*(4*q*q-aa*aa)*(bb*bb+4*q*q*amassp**2) c e3= print*,'e1,e1 ',e1*hc,e2*hc,epfermi*hc c epfermi=dsqrt(amassp**2+pfermi**2) efmin=max(e1,w+amassp-b,amassp,epfermi) efmax=w+dsqrt(amassp**2+pfermi**2)-b pmin=efmin pmax=efmax sekf=(pmax-pmin)/nstcos nsimps=(pmax-pmin)/sekf nsimps=2*(nsimps/2+1)+1 sekf=(pmax-pmin)/nsimps pmax=pmax+sekf/30.d0 c print*,'amin, amax, sa',pmin*hc,pmax*hc,sekf*hc if(pmin.gt.pmax) goto 12345 print*,'-------------------------------------------' print*,' at w=',w*hc,' MeV' write(*,*)' the 3-momentum transferred is q=',q*hc, ' MeV' print*,' the electron scattering angle is tet=',ang, ' degrees' print*,' the min. and max. value of Tf are:' print*,(pmin-amassp)*hc,(pmax-amassp)*hc suml0=0. sumt0=0. suml0t=0. sumt0t=0. sumtp0=0. sumtp0t=0. slwc0=0. slwc0t=0. stwc0=0. stwc0t=0. stpwc0=0. stpwc0t=0. wpc0=0. wpc0t=0. wpv0=0. wpv0t=0. do 1400 icharg=1,2 ffg=ffg0*z1num if(iabs(icharg).eq.2)ffg=ffg0*(z1mass-z1num) fac=4.d0 suml=0.d0 sumt=0.d0 sumtp=0.d0 sumlt=0.d0 sumtt=0.d0 sumtpt=0.d0 slwc=0.d0 slwct=0.d0 stwc=0.d0 stwct=0.d0 stpwc=0.d0 stpwct=0.d0 wpc=0.d0 wpct=0.d0 wpv=0.d0 wpvt=0.d0 kstart=1 do 1300 bet=pmin,pmax,sekf epf=bet epf2=epf**2 c print*,epf2-amassp**2 akpf=dsqrt(epf2-amassp**2) tpf=epf-amassp akpf2=akpf*akpf isimps=isimps+1 c ekf=eki-w Epi=-w+epf c printi,'epi+b-amassp ',epi+b-amassp,epi+b,epf,w akr=dsqrt(dabs((epi+b)**2-amassp*amassp)) akr=akr+1.d-11 pm=akr pm2=pm*pm c print*,'ekf ',ekf*hc c c calculo de w y Epi c w2=w*w Epi2=Epi*Epi c c calculo de wbar y ebar (energias iniciales si proton libre) c c c empiezan a calcularse las variables en el sistema en que pm=0 c ojo, todo calculado SIN RECOIL y en el sistema LABORATORIO c c Epibar=dsqrt(pm2+amassp*amassp) wbar=(Epf-Epibar) w2bar=wbar*wbar c c energia electron final c ekf2=ekf*ekf facos=2.0d0 sigacos=0.0d0 ichcos=0 sigacost=0.0d0 2393 format(/' e lept i=',f9.3,' e lept f=',f9.3,' akpf =',f9.3, , ' MeV'/ , ' Epf =',f9.3,' w =',f9.3,' pm =',f9.3, , ' MeV'/ , ' Epibar =',f9.3,' Tb =',f9.3,' Epi =',f9.3, , ' MeV'/ , ' Wbar =',f9.3,' eki =',f9.3,' ekf =',f9.3, , ' MeV'/ , ' Emiss =',f9.3,' q =',f9.1,' cosgam =',f9.6, , ' cosm =',f9.6 /) c c cosif=cosnu cosf=q2+akpf*akpf-pm*pm cosf=cosf/2.d0/q/akpf cosm=akpf/akr*cosf-q/akr if(imp.eq.1)then write(*,2393)eki*hc,ekf*hc,akpf*hc,epf*hc, ,w*hc,pm*hc,epibar*hc,tpf*hc,epi*hc,wbar*hc,eki*hc, , ekf*hc,emiss*hc,q*hc,cosf,cosm endif cosgam=cosf cos2gam=cosgam*cosgam singam=dsqrt(dabs(1.d0-cos2gam)) sin2gam=singam**2 c c los factores longitudinal/espacial etc c Epi2=Epi*Epi Epf2=Epf*Epf cosif=cosnu qnu2=w2-q2 qnu2bar=w2bar-q2 sin2if2=-qnu2/4.d0/eki/ekf cos2if2=1.d0-sin2if2 tg2if2=sin2if2/cos2if2 tgif2=dsqrt(dabs(tg2if2)) 818 format(11f7.1) c c rewrite with my notation... c some kinematical factors (only the terms that c survived to the azimuthal integration) c c Wl=qnu2*qnu2/q2/q2 Wpar=qnu2/q2/2.d0 * (-1.d0) +tg2if2 Wper=qnu2/q2/2.d0 * (-1.d0) +tg2if2 Wtt=0.0d0 Wlt=0.0d0 Wttp=-asig* 2.d0*tgif2 * dsqrt(dabs(tg2if2-qnu2/q2)) Wltp=0.0d0 c c some factors c free2=(Epf+Epibar)**2 c c se traen los factores de forma y la constante c de acoplo+propagador. Si ireac=0 es el c caso EM si es 1 es para neutrinos. Primero el caso EM c icharg es el indice del nucleon (1-proton, 2-neutron) c call ffactor(idip,icharg,0,qnu2,fff(1),eff,ggg,fcoup,f2s,g1s) C C MOTT'S CROSS SECTION in fm^-2 x units of fcoup (see in ffactor) C sigmott=4.d0*fcoup*ekf2*cos2if2 sm=sigmott if(iff.eq.2) fff(1)=0.d0 if(iff.eq.3) eff=0.d0 if(iff.eq.4) fff(1)=-eff if(iff.eq.5) then fff(1)=1.d-20 eff=1.d-20 ggg=1.d0 endif if(iff.eq.6)then ggg=0.d0 endif fff(2)=eff/2.d0/amassp f2=fff(2)**2 f1=fff(1)**2 f3=(eff+fff(1))**2 fff(3)=eff+fff(1) g1=ggg g12=g1*g1 gf=g1*(eff+fff(1)) pfz=akpf*cosgam pmz=pfz-q C write(57,'(7e11.3)')pm*hc,f1,f2,f3,g1,g12,gf c c f2 es (-kF2(q)/2Mp)**2, f1 es F1(q)**2, f3 es (f1+kF2)**2 c The free response functions (cc1 prescription) c rho2=free2*(f1-f2*qnu2bar)-q2*(f3 ) rho2=rho2+g12*(-4.d0*amassp**2 -q2+free2) djper2=-qnu2bar*f3 djper2=djper2+g12*(4.*amassp**2-qnu2bar) djpar2=akpf2*sin2gam*4.d0*(f1-f2*qnu2bar)+djper2 djpar2=djpar2+akpf2*sin2gam*4.d0*g12 djlt=0.d0 djtt=0.d0 djltp=0.0d0 djttp=-4.d0* gf * ((epf-epibar)* akpf * cosgam - epf* q) adjz2=(f1-f2*qnu2bar)*(pmz+pfz)**2- - f3*wbar**2+g12*(4*amassp**2-wbar**2+(pmz+pfz)**2) adjzr=g12*4*(epibar*pfz+epf*pmz)+ + f3*2.d0*(epibar-epf)*q+ + (f1-f2*qnu2bar)*2.d0*(epibar+epf)*(pmz+pfz) c c cambio de nombres... c rcc1=rho2 dltcc1=djlt percc1=djper2 parcc1=djpar2 djzcc1=djzz c c waw and the rest of structure functions have 1/fm^2 dimension c c c if icc<0 impose current conservation by removing J_z c if(icc.lt.0)then adjz2=w2/q2*rho2 adjzr=2.d0*w/q*rho2 endif rho2=rho2+w2/q2*adjz2-w/q*adjzr waw=( rho2+ wper*djper2 + wpar*djpar2 + wlt*djlt ) waw=waw+(wtt*djtt+wltp*djltp+wttp*djttp) sigcc1=waw*sigmott c c now the neutral current case c call ffactor(idip,icharg,1,qnu2,fwc(1),ewc,ggg,fcoupwc, , f2s,g1s) C if(iff.eq.2) fwc(1)=0.d0 if(iff.eq.3) ewc=0.d0 if(iff.eq.4) fwc(1)=-ewc if(iff.eq.5) then fwc(1)=1.d-20 ewc=1.d-20 ggg=1.d0 endif if(iff.eq.6)then ggg=0.d0 endif fwc(2)=ewc/2.d0/amassp f2wc=fwc(2)*fff(2) f1wc=fwc(1)*fff(1) f3wc=(ewc+fwc(1))*(eff+fff(1)) fwc(3)=ewc+fwc(1) g1=ggg g12=g1*g1 gf=g1*(eff+fff(1))/2.d0 C write(57,'(7e11.3)')pm*hc,f1,f2,f3,g1,g12,gf c c f2 es (-kF2(q)/2Mp)**2, f1 es F1(q)**2, f3 es (f1+kF2)**2 c The free response functions (cc1 prescription) c rho2wc=free2*(f1wc-f2wc*qnu2bar)-q2*(f3wc ) djper2wc=-qnu2bar*f3wc djpar2wc=akpf2*sin2gam*4.d0*(f1wc-f2wc*qnu2bar)+ + djper2wc djttpwc=-4.d0* gf * ((epf-epibar)* akpf * cosgam - epf* q) adjz2wc=(f1wc-f2wc*qnu2bar)*(pmz+pfz)**2- - f3wc*wbar**2 adjzrwc= + f3wc*2.d0*(epibar-epf)*q+ + (f1wc-f2wc*qnu2bar)*2.d0*(epibar+epf)*(pmz+pfz) c waw and the rest of structure functions have 1/fm^2 dimension c c c if icc<0 impose current conservation by removing J_z c if(icc.lt.0)then adjz2wc=w2/q2*rho2wc adjzrwc=2.d0*w/q*rho2wc endif rho2wc=rho2wc+w2/q2*adjz2wc-w/q*adjzrwc wawwc=rho2wc+ wper*djper2wc + wpar*djpar2wc wawwc=-(CA*wawwc-CV*wttp*djttpwc) sigwc=2.d0*wawwc*sigmott/dsqrt(fcoup)*dsqrt(fcoupwc) c c so sigcc1 has the units of fcoup x fm^(-4) c and in principle that means it is adimensional c c acumulate response functions c EM responses suml=suml+rho2*fac sumt=sumt+(djper2*wper+djpar2*wpar)*fac sumtp=sumtp+wttp*djttp*fac sumlt=sumlt+rho2 sumtt=sumtt+djper2*wper+djpar2*wpar sumtp=wttp*djttp c PV responses slwc=slwc+rho2wc*fac stwc=stwc+(djper2wc*wper+djpar2wc*wpar)*fac stpwc=stpwc+wttp*djttpwc*fac slwct=slwct+rho2wc stwct=stwct+djper2wc*wper+djpar2wc*wpar stpwct=stpwct+wttp*djttpwc c acumulate cross-sections wpc=wpc+sigcc1*fac wpv=wpv+sigwc*fac wpct=wpct+sigcc1 wpvt=wpvt+sigwc fac=6.d0-fac write(36,'(2f9.3,3e13.5,f11.5)')w*hc,pm*hc, , wpc,wpv c c keep values at starting point c if (kstart.eq.1)then a0l=rho2 a0t=wper*(djper2+djpar2) a0tp=wttp*djttp a0lwc=rho2wc a0twc=wper*(djper2wc+djpar2wc) a0tpwc=wttp*djttpwc w0l=sigcc1 w0t=sigwc kstart=0 endif 1300 continue c c include the Fermi gas normalization and the factor c 1/q, and along with sekf (in fm) that makes the c cross section in fm^3 c factor=ffg*(2.d0*pi)/q/4.d0 factor=factor/hc c c ffg is the corresponding normalization for RFG c factor is in fm^4 so esig is in the units of fcoup c c c Simpsons 3/4 integrals c suml=suml0+factor*(suml-rho2-a0l)*sekf/3.d0 sumt=sumt0+factor*(sumt-wper*(djper2+djpar2)-a0t)*sekf/3.d0/wper sumtp=sumtp0+factor*(sumtp-wttp*djttp-a0tp)*sekf/3.d0/wttp slwc=slwc0+factor*(slwc-rho2wc-a0lwc)*sekf/3.d0 stwc=stwc0+factor*(stwc-wper*(djper2wc+djpar2wc)-a0twc)* 1 sekf/3.d0/wper stpwc=stpwc0+factor*(stpwc-wttp*djttpwc-a0tpwc)*sekf/3.d0/wttp wpc=wpc0+factor*(wpc-sigcc1-w0l)*sekf/3.d0 wpv=wpv0+factor*(wpv-sigwc-w0t)*sekf/3.d0 c c same as before but now with trapezoidal's rule c sumlt=suml0t+factor*(sumlt-0.5*rho2-0.5*a0l)*sekf sumtt=sumt0t+factor*(sumtt-0.5*wper*(djper2+djpar2)- 1 0.5*a0t)*sekf/wper sumtpt=sumtp0t+factor*(sumtpt-0.5*wttp*djttp-0.5*a0tp) 1 *sekf/wttp slwct=slwc0t+factor*(slwct-0.5*rho2wc-0.5*a0lwc)*sekf stwct=stwc0t+factor*(stwct-0.5*wper*(djper2wc+djpar2wc)- 1 0.5*a0twc)*sekf/wper stpwct=stpwc0t+factor*(stpwct-0.5*wttp*djttpwc- 1 0.5*a0tpwc)*sekf/wttp wpct=wpc0t+factor*(wpct-0.5*sigcc1-0.5*w0l)*sekf wpvt=wpv0t+factor*(wpvt-0.5*sigwc-0.5*w0t)*sekf if(icharg.eq.1)then write(37,'(f7.3,6e12.4)')w*hc,wpv/wpc,wpvt/wpct endif if(icharg.eq.2)then write(38,'(f7.3,6e12.4)')w*hc,1./((wpc-wpc0)/(wpv-wpvt)), , 1./( (wpct-wpc0t)/(wpvt-wpv0t)) print*,'and the assymetry is',wpv/wpc,ang c c now printing the asymmetry and cross-section in fm^3 c c asymmetry calculation WPCTOT=suml+wper*sumt WPVTOT=-(CA*(slwc+wper*stwc)-CV*wttp*stpwc) ASYMM=2.d0*dsqrt(fcoupwc/fcoup)*(WPVTOT/WPCTOT) c WPCTOTt=sumlt+wper*sumtt WPVTOTt=-(CA*(slwct+wper*stwct)-CV*wttp*stpwct) ASYMMt=2.d0*dsqrt(fcoupwc/fcoup)*(WPVTOTt/WPCTOTt) c c c c c write(39,'(f7.1,6e12.4)')w*hc,wpv/wpc,wpvt/wpct,asymm,asymmt c c write(40,'(f9.2,e13.5,2e10.3,e13.5,2e10.3)')q*hc, ,wpv/wpc,wpv,wpc,wpvt/wpct,wpvt,wpct write(50,'(f7.1,6e12.4)')w*hc,suml/wl,sumlt/wl, 1 sumt,sumtt,sumtp,sumtpt write(51,'(f7.1,6e12.4)')w*hc,ca*slwc/wl,ca*slwct/wl, 1 ca*stwc,ca*stwct,2*cv*stpwc,2*cv*stpwct endif suml0=suml suml0t=sumlt sumt0=sumt sumt0t=sumtt sumtp0=sumtp sumtp0t=sumtpt slwc0=slwc slwc0t=slwct stwc0=stwc stwc0t=stwct stpwc0=stpwc stpwc0t=stpwct wpc0=wpc wpc0t=wpct wpv0=wpv wpv0t=wpvt 1400 continue 1 continue 12345 continue close(39) close(37) close(38) close(46) close(47) close(13) close(57) close(67) close(68) end subroutine ffactor(idip,icharg,ireac,qnu2,ef1,ef2,eg1,fcoup, , f2s,g1s) implicit real*8 (a-h,o-z) real*8 eff1(2),eff2(2), , ff1(2),ff2(2),f1(0:8),f2(0:8),g1(0:8),gg1(2) akp=1.79285d0 akn=-1.91304d0 alph=1.d0/137.04 pi=3.1415926 c c corregir el valor aqui tambien! c sin2w=0.227 c c Gf in MeV^(-2) c Gf=1.1664e-11 sqr3=sqrt(3.0) amass=939.d0 hc=197.3d0 am2v=840.d0**2 am2a=1030.d0**2 c f2s=-0.21 c g1s=-0.19 c f2s=0. c g1s=0. f1s=0. f1(0)=1.+f1s f2(0)=f2s+(akp+akn) g1(0)=g1s+0.55/3. f1(3)=0.5 f2(3)=0.5*(akp-akn) g1(3)=0.63 f1(8)=0.5*sqr3 f2(8)=0.5*sqr3*(akp+akn) g1(8)=0.275/sqr3 z0v=-0.5 z0a=0.5 z3v=1.-2*sin2w z3a=-1. z8v=1./sqr3*(1.-2*sin2w) z8a=-1./sqr3 c print*,z0v,z3v,z8v,f2(0),f2(3),f2(8) c write(10,*)'proton-nu proton anu neutron nu neutron-anu' Q2=-qnu2*hc*hc dip=1.+q2/am2v dip=1./dip/dip dips=1.+Q2/am2a dips=1./dips/dips tau=Q2/4./amass/amass c------------------------------------ c c factores de forma EM del neutron y proton c if (icharg.eq.1)then anom=akp ge0=1.0 else anom=akn ge0=0.0 endif Gm0=ge0+anom if(idip.eq.0)then dip=1.d0 dips=1.d0 endif gm=gm0*dip ge=ge0*dip if(icharg.eq.2)then ge=-anom*tau/(1.+5.6*tau)*dip endif eff2(icharg)=(gm-ge)/(1.+tau)/anom eff1(icharg)=(ge+tau*gm)/(1.+tau) gmt=eff1(icharg)+anom*eff2(icharg) get=eff1(icharg)-tau*anom*eff2(icharg) f20=f2(0)/(1.+tau)*dip f23=f2(3)/(1.+tau)*dip f28=f2(8)/(1.+tau)*dip f10=f1(0)*dip+tau*f20 f13=f1(3)*dip+tau*f23 f18=f1(8)*dip+tau*f28 g10=g1(0)*dips g13=g1(3)*dips g18=g1(8)*dips c c fac is for weak case (Gf/sqrt(2))^2 so for the moment is c in MeV^(-4) c fac=(gf/sqrt(2.d0)/4.d0/pi)**2 ff1(1)=z0v*f10+z3v*f13+z8v*f18 ff1(2)=z0v*f10-z3v*f13+z8v*f18 ff2(1)=z0v*f20+z3v*f23+z8v*f28 ff2(2)=z0v*f20-z3v*f23+z8v*f28 gg1(1)=z0a*g10+z3a*g13+z8a*g18 gg1(2)=z0a*g10-z3a*g13+z8a*g18 c print*,'ff1(1),ff1(2),gg1(1),gg1(2)',ff1(1),ff1(2),gg1(1),gg1(2) c print*,'ff2 ',ff2(1),ff2(2) iii=icharg c write(13,'(i4,5e13.5)')iii,q2,ff1(iii),ff2(iii),gg1(iii) c i=1, proton c i=2, neutron if(ireac.ne.0)then ef1=ff1(icharg) ef2=ff2(icharg) eg1=gg1(icharg) c c neutral case, fac is in MeV-4 convert to fm^4 c with a factor c (hc)^4 c fcoup=fac*(hc)**4 else eg1=0.d0 ef1=eff1(icharg) ef2=eff2(icharg)*anom c c em case: fcoup is in fm^4 c fcoup=alph*alph/qnu2/qnu2 endif write(57,'(4f10.5,f10.6,2f10.5)')qnu2,ef1,ef2,eg1,fcoup, , dip,gmt/get end