	PROGRAM PHITOJ
c
c  R.T. AZUAH		JULY 93
c ######################################################################
c Converts a tof scan from a single detector which is in S(phi,w) to 
c J(Y). Assumes the IA holds in which case any tof scan which cuts through 
c the recoil line should give us J(Y). Note that J(Y) (not J(Q,Y) ) does 
c not depend on Q and w separately but on both of them through Y.
c
c ** Reads SMB SPEAL output (note shift in phi bins here)
c
	integer*4 n_phi,n_e,idec,in_len,strim,i,n_y,j,maxsplit,rebin_data,n_q
	real*4 phi(300),ene(1000),sjy(1000),ejy(1000),sums(1000),sume(1000)
	real*4 sqe(1000),eqe(1000),M,Ei,pi,ki,kf,const1,const2,const3,ang(300)
	real*4 y_o(1000),Q(1000),ymin,ymax,d2r,r2d,const4,ybin(1000)
	real*8 dy,y(1000)
	real*4 sqe2d(300,1000),eqe2d(300,1000),sjy2d(300,1000),ejy2d(300,1000)
	real*4 Qmin,Qmax,dQ,phibin(300),Qb,E
	character*50 filein,fileout,ext*4
	character*50 text1,text2,input*10
	write (*,'(T5,A,$)') 'Give Input File name:: '
	read (*,'(A)') filein
	write (*,'(T5,A,$)') 'Give Ei and rel at mass:: '
	read (*,*) Ei,M
	write (*,*) 'Give min, max and interval of Y :: '
	read (*,*) ymin,ymax,dy
	write (*,*) 'Give min, max and interval of Q :: '
	read (*,*) Qmin,Qmax,dQ
	in_len = strim(filein)
	fileout=filein(1:in_len-3)//'ijy'
	open (unit=10,file=filein,status='old')
	const1=2.07219
	const2=const1*1.008665/M
	const3=const1*1.008665*2.0
	pi=4.0*ATAN(1.0); d2r = pi/180.0; r2d = 1/d2r
	ki=SQRT(Ei/const1)
c
c Read phi and energy grids from input file
c Write phi (>=120 deg) and Y grid to output file
c
	read (10,*) n_phi,n_e
	read (10,'(A)') text1
	read (10,'(8G10.4)') (phi(i),i=1,n_phi)
	read (10,'(A)') text2
	read (10,'(8G10.4)') (ene(i),i=1,n_e+1)
	n_y = 1
	y(n_y) = ymin
	ybin(n_y) = y(n_y) - 0.5*dy 
	DO WHILE ( y(n_y) .LT. (ymax-0.5*dy) )
	   n_y = n_y + 1
	   y(n_y) = y(n_y - 1) + dy
	   ybin(n_y) = ybin(n_y -1) + dy
	END DO
	ybin(n_y + 1) = ybin(n_y) + dy
	n_y = n_y - 1
c Create required Q-Grid and the corresponding Phi-Grid calculated at the recoil Energy for each Q
c Note that the calculated Phibin are histogrm values, not the true values.
	n_q = 1
	Q(n_q) = Qmin
	Qb = Qmin - 0.5*dQ
	E = const2*Qb**2
	phibin(n_q) = 2*ki**2 - E/const1 - Qb**2
	phibin(n_q) = phibin(n_q) / (2*ki*SQRT(ki*ki - E/const1))
	phibin(n_q) = r2d*ACOS(phibin(n_q))
	DO WHILE ( Q(n_q) .LT. (Qmax-0.5*dQ) )
	   n_q = n_q + 1
	   Q(n_q) = Q(n_q - 1) + dQ
	   Qb = Qb + dQ
	   E = const2*Qb**2
	   phibin(n_q) = 2*ki**2 - E/const1 - Qb**2
	   phibin(n_q) = phibin(n_q) / (2*ki*SQRT(ki*ki - E/const1))
	   phibin(n_q) = r2d*ACOS(phibin(n_q))
	END DO
	Qb = Qb + dQ
	E = const2*Qb**2
	phibin(n_q+1) = 2*ki**2 - E/const1 - Qb**2
	phibin(n_q+1) = phibin(n_q+1) / (2*ki*SQRT(ki*ki - E/const1))
	phibin(n_q+1) = r2d*ACOS(phibin(n_q+1))
c Open Output file and write the Q and E headers
	open (unit=20,file=fileout,status='unknown')
	write (20,'(2I5,G10.4)') n_q,n_y+1,ei
	write (20,'(A)') '### Q Grid'
	write (20,'(8G10.4)') (Q(i),i=1,n_q)
	write (20,'(A)') '### Y Grid'
	write (20,'(8G10.4)') (y(i),i=1,n_y+1)
c	write (*,'(8G10.4)') (y(i),i=1,n_y+1)
c	write (*,'(8G10.4)') (ybin(i),i=1,n_y+1)

c
c Calculate the histogram values of the input angles.
c
	DO i=1,n_phi-1
	   ang(i) = 0.5*(3*phi(i) - phi(i+1))
c	   write (*,*) ' Proper angle = ',ang(i)
	END DO
	ang(n_phi) = 0.5*(phi(n_phi-1)+phi(n_phi))
c	write (*,*) ' Proper angle = ',ang(n_phi)
	ang(n_phi+1) = 0.5*(3*phi(n_phi)-phi(n_phi-1))
c	write (*,*) ' Proper angle = ',ang(n_phi+1)
c
c Loop through blocks, transform data to J(Y) and rebin to requested Y-Grid
c
	DO i=1,n_phi
	   read (10,'(A)') text1
	   read (10,'(8G10.4)') (sqe2d(i,j),j=1,n_e) 
	   read (10,'(A)') text2
	   read (10,'(8G10.4)') (eqe2d(i,j),j=1,n_e)
	   DO j=1,n_e+1
	      kf=SQRT((Ei-ene(j))/const1)
	      Q(j) = SQRT(ki**2 + kf**2 - 2*ki*kf*COS(phi(i)*d2r))
	      const4 = Q(j)*const3/M
	      y_o(j) = (ene(j)-const2*Q(j)**2)/const4
	      sqe(j) = sqe2d(i,j)*const4
	      eqe(j) = (eqe2d(i,j)*const4)**2
	   END DO
	   IF (rebin_data(n_e,y_o,sqe,eqe,n_y,ybin,sjy,ejy,maxsplit) .EQ. 0) THEN
	      write (*,*) ' ERROR: REBIN_DATA Routine failed at n_phi = ',i
	      STOP
	   END IF
	   DO j=1,n_y
	      sjy2d(i,j) = sjy(j)
	      ejy2d(i,j) = ejy(j)
	   END DO
c	   write (20,'(A)') '### J(Y)'
c	   write (20,'(8G10.4)') (sjy(j),j=1,n_y)
c	   write (20,'(A)') '### Errors'
c	   write (20,'(8G10.4)') (SQRT(ejy(j)),j=1,n_y)
	END DO
c
c Rebin data to requested Q-Grid and output results to output file
c

	DO j = 1,n_y
	   DO i = 1,n_phi
	      sqe(i) = sjy2d(i,j)
	      eqe(i) = ejy2d(i,j)
	   END DO
	   IF (rebin_data(n_phi,ang,sqe,eqe,n_q,phibin,sjy,ejy,maxsplit) .EQ. 0) THEN
	      write (*,*) ' ERROR: REBIN_DATA Routine failed at n_y = ',j
	      STOP
	   END IF
	   DO i=1,n_q
	      sjy2d(i,j) = sjy(i)
	      ejy2d(i,j) = ejy(i)
	   END DO
	END DO
c
c Write Output file
c


	DO i = 1,n_q
	   write (20,'(A)') '### J(Y)'
	   write (20,'(8G10.4)') (sjy2d(i,j),j=1,n_y)
	   write (20,'(A)') '### Errors'
	   write (20,'(8G10.4)') (SQRT(ejy2d(i,j)),j=1,n_y)
	END DO

	STOP
	END














