! qvec.f
! 
! gh, 07.07.31
!
! program to calculate omega and q values for given Q2, Pmiss, following
! Steffen's idea given below.
!
! In parallel kinematics we have for the two-body breakup:
!
! omega + M(4He)  = sqrt(p2 + m2) + sqrt(pr2 + M(3H)3)
! q = p + pr
!
! where p is the proton momentum and pr is the missing momentum; omega
! and q are the energy and momentum transfers. From here we get:
!
! omega + M(4He)  = sqrt((q - pr)2 + m2) + sqrt(pr2 + M(3H)3)

! Thus with omega and q fixed (electron arm fixed) the missing momentum
! is fixed, too. If we only change the momentum of the proton arm we
! just probe different (omega, q) values which are allowed by the
! acceptance of the electron spectrometer. We should in my opinion rather
! use the Q2 to fix the omega and q relation in the above equation. Then,
! the above equation tells us what omega and q needs to be if we want to
! reach a certain pr. We get the proton momentum, too, and should set the
! hadron arm central momentum to that value.

      implicit none

      real*8 xmp,xmn,xmd,xmt,xm4he,xm3he,mtarget,mrecoil
      real*8 qsq(7),ebeam(7),pmiss(5),omega
      real*8 omega1,omega2,q,qnew,qold,err
      real*8 pshms,stheta2,thshms,phms,tanhms,thhms
      real*8 q2chk1,q2chk2
      integer i,j

c stuff for linux_suppl.inc
      real*8 sind, cosd, tand
      real*8 asind, acosd, atand, atan2d
      external sind, cosd, tand
      external asind, acosd, atand, atan2d

! Note Q2=2.47 is where the 2-gamma expt will measure Gep/Gmp to <1%.
      data qsq  /1.0,  1.15, 1.3,  2.0,  2.1,  2.2,  2.47/
      data ebeam/2.25, 2.25, 2.25, 4.40, 4.40, 4.40, 4.40/
      data pmiss/0.250,0.140,0.00,-0.140,-0.250/

      xmp   = 0.93827
      xmn   = 0.93957
      xmd   = 1.87561
      xmt   = 2.80892
      xm3he = 2.80842
      xm4he = 3.72742

      mtarget=xmd
      mrecoil=xmn
*      mtarget=xm4he
*      mrecoil=xm3he

! the equations are simple in concept but involve many square roots.
! so i solve them iteratively

      do i=1,7    ! Q2 loop
 5       continue
         write(6,*)
         write(6,10)qsq(i),ebeam(i)
 10      format('  **** Q2=',f7.3,'   Ebeam=',f5.3)

         do j=1,5 ! PM loop

            if (qsq(i).gt.1.3 .and. abs(pmiss(j)).gt.0.150) goto 200
            if (qsq(i).gt.1.6 .and. abs(pmiss(j)).gt.0.050) goto 200

            q=1.0
            err=0.1

 50         omega1=q**2-qsq(i)
            omega2=( sqrt((q-pmiss(j))**2+xmp**2) 
     1           +sqrt(pmiss(j)**2+mrecoil**2) - mtarget )**2

            if (abs(omega1-omega2).lt.err .and.err.lt.1.e-6) goto 100

            if (abs(omega1-omega2).lt.err) err=err/2.
            
            if (omega1.gt.omega2) then   ! q too large
               qnew=q-err
            else
               qnew=q+err
            endif
c           write(6,*)q,omega1,omega2,qnew

            qold=q
            q=qnew
            goto 50

 100        omega=sqrt(omega1)

            q2chk1=q**2-omega**2

            pshms=ebeam(i)-omega
            stheta2=qsq(i)/(4.*ebeam(i)*pshms)
            thshms=2.*asind(sqrt(stheta2))

            q2chk2=4.*ebeam(i)*pshms*(sind(thshms/2.))**2

c           write(6,*)q2chk1,q2chk2

            phms=q-pmiss(j)
            tanhms=(pshms*sind(thshms))/(ebeam(i)-pshms*cosd(thshms))
            thhms=atand(tanhms)

            if (thshms.gt.40.0) then
               write(6,1001)thshms,pmiss(j)
 1001          format(' ThSHMS too large ',f6.2,f8.3)
               ebeam(i)=ebeam(i)+0.050
               goto 5
            endif

            write(6,101)pmiss(j),omega,q,pshms,thshms,phms,thhms
 101        format(' PM=',f7.3,'    Omega=',f5.3,' q=',f5.3,
     1           '   PSHMS=',f5.3,' ThSHMS=',f6.2,
     2           '   PHMS=',f5.3,' ThHMS=',f6.2)

 200     enddo
      enddo
            
      end

      include 'linux_suppl.inc'