program hadji c main program to run routine ffhadj c Hadjimichael c The output if file "rhop.dat" c The contents depends of nopt below c nopt=1 :proton and neutron electric ff c 2 :He3 and H3 charge ff c 3 :He3 and H3 magnetic ff c 4 :Charge He3 form factor Implicit Real*8 (a-h,o-z) q2mev=0.0 Open(6,FILE='rhop.dat', Status='Unknown') nif=1 !it is nifty(19) no MEC nopt=4 q=0.05 Do i=1,200 q2mev=(q*197.329)**2 !to give (Mev/c)**2 units qfer2=q*q If(nopt.eq.1) Then Call ffpn(qfer2,fp,fn) Write(6,*)q*q,fp,fn q=q+0.05 End If If((nopt.eq.2).or.(nopt.eq.3))Then Call ffmec(qfer2,fche3,fch3,fmhe3,fmh3,nif) if(nopt.eq.2)Write(6,*) q,abs(fche3),abs(fch3) if(nopt.eq.3)Write(6,*) q,abs(fmhe3),abs(fmh3) q=q+0.05 End If if(nopt.eq.4) Then Call ffhadj(q2mev,rp,rn,rpsp,rnsp,nif) Write(6,100) qfer2,abs(rp),abs(rn),abs(rpsp),abs(rnsp) q=q+0.024 End If End Do 100 Format(1h , 5(e13.5)) Close(6) Stop 'finished' End Subroutine ffhadj(q2mev,rp,rn,rpsp,rnsp,nif) c************************************************************************** c * c calculates the matter and spin form factors of the protons and neutron * c in helium-3 from the charge and magnetic form factors of he-3 and h-3 * c as calculate by hadjimichael et. al. in the impulse approximation. * c * c thus no meson currents are in the charge and magnetic form factors. * c * c************************************************************************** Implicit Real*8 (a-h,o-z) data up,un/2.79278,-1.91315/ c************************************************************************** c * c q2mev must be converted to inverse Fermis squared * c * c************************************************************************** qfer2 = q2mev/(197.329**2) c c c************************************************************************** c * c now calculate the form factor of the proton * c * c************************************************************************** c c Call ffpn(qfer2,fp,fn) c c c************************************************************************** c * c now calculate the charge and magnetic form factors of he-3 and h-3 * c * c************************************************************************** c c Call ffmec(qfer2,fche3,fch3,fmhe3,fmh3,nif) c c c************************************************************************** c * c now calculate matter and spin form factors using formulas of mach that * c relate these to the charge and magnetic form factors of he-3 and h-3. * c LPOTT eqs. 80- 83 * c************************************************************************** c c rp = fche3/fp rn = fch3/fp denom = fp*((up**2)-(un**2)) rpsp = ((up*un)*(fmhe3-fmh3))/(2*denom) up2 = up**2 un2 = un**2 rnsp = (up2*fmh3-un2*fmhe3)/denom Return End Subroutine ffmec(q2f2,ff1,ff2,ff3,ff4,n) c************************************************************************** c * c this Subroutine calculates the charge and magnetic form factors of * c he-3 and h-3 using the impulse approximation of hadjimichael et. al. * c (phys rev c vol 27, no.2, p831, feb 1983). * c * c************************************************************************** c****************************************************************************** c * c up dated on sept 9, 1984 to enable the calculation of the form factors * c from hadjimicheal et al. using the total charge, and magnetic form factors * c * c****************************************************************************** Implicit Real*8 (a-h,o-z) Real*8 lpt Dimension formfs(4,20),y(4),q2hadj(20) Dimension form1(20),form2(20),form3(20),form4(20) Dimension ftot1(20),ftot2(20),ftot3(20),ftot4(20) Dimension fhadj(4),a2(4) Dimension a2tot(4) c c data q2hadj/ 0.0, 0.25, 1.0, 2.25, 4.0, 6.25, 9.0, 12.25, 16.0, 1 20.25, 25.0, 30.25, 36.0, 42.25, 49.0, 56.25, 2 64.0, 72.25, 81.0, 90.25/ c data a2/ 39.474, 3.012, 16.667, 18.204/ c data a2tot/28.667, 21.928, 26.864, 83.324/ c data form1/1.0,0.865,0.571,0.308,0.141,0.552d-01,0.174d-01, 10.332d-02,-0.919d-03,-0.150d-02,-0.112d-02,-0.673d-03,-0.349d-03, 2-0.162d-03,-0.659d-04,-0.223d-4,-0.445d-05,0.133d-05,0.251d-05, 30.203d-5/ c data form2/1.0,0.868,0.604,0.347,0.170,0.719d-1,0.252d-01, a 0.612d-02, 1-0.294d-03,-0.171d-02,-0.152d-02,-0.102d-02,-0.595d-03, 2-0.316d-03,-0.154d-03,-0.696d-04,-0.285d-04,-0.101d-04,-0.261d-05, 3-0.691d-06/ c data form3/1.0,0.843,0.522,0.246,0.873d-01,0.163d-01,-0.719d-02, 1-0.107d-01,-0.818d-02,-0.494d-02,-0.257d-02,-0.116d-02,-0.419d-03, 2-0.826d-04,0.480d-04,0.798d-04,0.735d-04,0.554d-04,0.379d-04, 30.241d-04/ c data form4/1.0,0.861,0.566,0.298,0.130,0.462d-01,0.113d-01, a-0.313d-03, 1-0.274d-02,-0.232d-02,-0.143d-02,-0.740d-03,-0.321d-03,-0.107d-03, 2-0.104d-04,0.233d-04,0.296d-04,0.253d-04,0.188d-04,0.125d-04/ c data ftot1/1.0,0.849,0.556,0.294,0.129,0.445d-01,0.890d-02, a-0.276d-02,-0.480d-02,-0.376d-02,-0.232d-02,-0.123d-02, b-0.549d-03,-0.185d-03,-0.124d-04,0.532d-04,0.679d-04, c0.608d-04,0.475d-04,0.344d-04/ c data ftot2/1.0,0.868,0.604,0.347,0.169,0.705d-01,0.236d-01, a0.462d-02,-0.162d-02,-0.279d-02,-0.235d-02,-0.163d-02, b-0.102d-02,-0.611d-03,-0.355d-03,-0.206d-03,-0.121d-03, c-0.735d-04,-0.462d-04,-0.303d-04/ c data ftot3/1.0,0.863,0.576,0.317,0.151,0.625d-01, a0.221d-01,0.587d-02,0.358d-03,-0.990d-03,-0.975d-03, b-0.665d-03,-0.376d-03,-0.183d-03,-0.711d-04,-0.153d-04, c0.840d-05,0.148d-04,0.138d-04,0.978d-05/ c data ftot4/1.0,0.874,0.599,0.345,0.172,0.770d-01, a0.307d-01,0.108d-01,0.301d-02,0.349d-03,-0.343d-03, b-0.401d-03,-0.287d-03,-0.171d-03,-0.882d-04,-0.390d-04, c-0.132d-04,-0.160d-05,0.247d-05,0.276d-05/ c c c************************************************************************** c * c If q2f2 outside of table then extrapolate with * c * c f(q2) = f(q02)exp -(q2-q02)/a2 * c * c where f(q02) = last point in hadjimichaels table (ie. q02=90.25 or * c q = 9.5) * c * c************************************************************************** c c If (n .eq. 2) GoTo 10 Do 70 j=1,20 formfs(1,j) = form1(j) formfs(2,j) = form2(j) formfs(3,j) = form3(j) formfs(4,j) = form4(j) 70 Continue GoTo 15 10 Continue Do 75 j=1,20 formfs(1,j) = ftot1(j) formfs(2,j) = ftot2(j) formfs(3,j) = ftot3(j) formfs(4,j) = ftot4(j) 75 Continue 15 Continue If (q2f2 .lt. 90.25) GoTo 40 lpt = 90.25 If (n .eq.2) GoTo 25 Do 20 i=1,4 q2a2 = (q2f2-lpt)/a2(i) If (q2a2 .gt. 150.) q2a2 = 150. fhadj(i) = formfs(i,20)*exp(-q2a2) 20 Continue GoTo 30 25 Continue Do 26 i=1,4 q2a2 =(q2f2-lpt)/a2tot(i) If (q2a2 .gt. 150.) q2a2 = 150. fhadj(i) = formfs(i,20)*exp(-q2a2) 26 Continue 30 Continue ff1 = fhadj(1) ff2 = fhadj(2) ff3 = fhadj(3) ff4 = fhadj(4) Return 40 Continue If ((q2f2 .gt. 81.0) .or. (q2f2 .lt. 0.25)) GoTo 50 c c c************************************************************************* c * c for q2f2 less than 81.0 use four point lagrangian interpolation. * c * c************************************************************************* c c Call lagrng(q2f2,q2hadj,y,formfs,20,4,4,20,4) ff1 = y(1) ff2 = y(2) ff3 = y(3) ff4 = y(4) Return c c 50 Continue c c c************************************************************************** c * c for q2f2 less than 90.25 but greater than 81.0 must use 2-point * c lagrangian interpolation since in between the last two data points * c of table. * c * c************************************************************************** c c Call lagrng(q2f2,q2hadj,y,formfs,20,4,2,20,4) ff1 = y(1) ff2 = y(2) ff3 = y(3) ff4 = y(4) Return End Subroutine ffpn (q2,fcp,fcn) c *** see RH Landau, Program lpott, 1981 c proton charge ff (also p mag and n mag If nrom(0)=1) c neutron charge form factor c q in f-1 Implicit Real*8 (a-h,m,o-z) c data x,un,mn/18.2,-1.91315,4.7548/ fcp = 1./((1.+q2/x)**2) fcn = -(un*q2*0.25/mn/mn)*fcp Return End Subroutine lagrng (x,arg,y,val,ndim,nfs,nptx,maxarg,maxfs) c *** see RH Landau, Program lpott, 1981 c lagrange interpolation,unequally spaced points c npts=2,3,4,5,6, nfs functions(y*s) simultaneously interpltd c x= value of argument, arg is the tabulated x*s c y= a vector of interpolad functions, from tabulted val c ndim= Dimension of table c nfs= = of functions simult interpolated c maxarg and maxfs are maximum values of subscripts,ie Dimensions Implicit Real*8 (a-h,o-z) Dimension arg(maxarg), val(maxfs,maxarg), y(maxfs) c -----find x0 the Closest point to x npts = nptx c 2=st ni = 1 nf = ndim 10 If ((x.le.arg(ni)).or.(x.ge.arg(nf))) GoTo 30 If ((nf-ni+1).eq.2) GoTo 70 nmid = (nf+ni)/2 If (x.gt.arg(nmid)) GoTo 20 nf = nmid GoTo 10 20 ni = nmid GoTo 10 c ------ x is one of the tabltd values 30 If (x.le.arg(ni)) GoTo 60 nn = nf 40 nused = 0 Do 50 n=1,nfs y(n) = val(n,nn) 50 Continue Return 60 nn = ni GoTo 40 c ------- 2 pts left, choose smaller one 70 n0 = ni nn = npts-1 GoTo (140,110,100,90,80), nn 80 Continue If (((n0+3).gt.ndim).or.((n0-2).lt.1)) GoTo 90 nused = 6 GoTo 150 90 Continue If ((n0+2).gt.ndim) GoTo 110 If ((n0-2).lt.1) GoTo 100 nused = 5 GoTo 150 100 Continue If (((n0+2).gt.ndim).or.((n0-1).lt.1)) GoTo 110 nused = 4 GoTo 150 110 If ((n0+1).lt.ndim) GoTo 130 c n0=ndim,special case 120 nn = ndim GoTo 40 130 nused = 3 If ((n0-1).lt.1) nused = 2 GoTo 150 140 If ((n0+1).gt.ndim) GoTo 120 nused = 2 150 Continue c at least 2 pts left y0 = x-arg(n0) y1 = x-arg(n0+1) y01 = y1-y0 c0 = y1/y01 c1 = -y0/y01 If (nused.eq.2) GoTo 160 c at least 3 pts ym1 = x-arg(n0-1) y0m1 = ym1-y0 ym11 = y1-ym1 cm1 = -y0*y1/y0m1/ym11 c0 = c0*ym1/y0m1 c1 = -c1*ym1/ym11 If (nused.eq.3) GoTo 180 c ------at least 4 pts y2 = x-arg(n0+2) ym12 = y2-ym1 y02 = y2-y0 y12 = y2-y1 cm1 = cm1*y2/ym12 c0 = c0*y2/y02 c1 = c1*y2/y12 c2 = -ym1*y0*y1/ym12/y02/y12 If (nused.eq.4) GoTo 200 c at least 5 pts ym2 = x-arg(n0-2) ym2m1 = ym1-ym2 ym20 = y0-ym2 ym21 = y1-ym2 ym22 = y2-ym2 cm2 = ym1*y0*y1*y2/ym2m1/ym20/ym21/ym22 cm1 = -cm1*ym2/ym2m1 c0 = -c0*ym2/ym20 c1 = -c1*ym2/ym21 c2 = -c2*ym2/ym22 If (nused.eq.5) GoTo 220 c at least 6 pts y3 = x-arg(n0+3) ym23 = y3-ym2 ym13 = y3-ym1 y03 = y3-y0 y13 = y3-y1 y23 = y3-y2 cm2 = cm2*y3/ym23 cm1 = cm1*y3/ym13 c0 = c0*y3/y03 c1 = c1*y3/y13 c2 = c2*y3/y23 c3 = ym2*ym1*y0*y1*y2/ym23/ym13/y03/y13/y23 GoTo 240 160 Continue Do 170 n=1,nfs y(n) = c0*val(n,n0)+c1*val(n,n0+1) 170 Continue GoTo 260 180 Continue Do 190 n=1,nfs y(n) = cm1*val(n,n0-1)+c0*val(n,n0)+c1*val(n,n0+1) 190 Continue GoTo 260 200 Continue Do 210 n=1,nfs y(n) = cm1*val(n,n0-1)+c0*val(n,n0)+c1*val(n,n0+1)+c2*val(n,n0+ 1 2) 210 Continue GoTo 260 220 Continue Do 230 n=1,nfs y(n) = cm2*val(n,n0-2)+cm1*val(n,n0-1)+c0*val(n,n0)+c1*val(n,n0 1 +1)+c2*val(n,n0+2) 230 Continue GoTo 260 240 Continue Do 250 n=1,nfs y(n) = cm2*val(n,n0-2)+cm1*val(n,n0-1)+c0*val(n,n0)+c1*val(n,n0 1 +1)+c2*val(n,n0+2)+c3*val(n,n0+3) 250 Continue 260 Return End