1 gaskelld 1.1 	subroutine mc_hrsr_recon (delta_p,delta_t,delta_phi,y_tgt,fry)
  2              
  3              C+______________________________________________________________________________
  4              !
  5              ! MC_HRSR_RECON : Reconstruct target quantities from tracks.
  6              !		   This subroutine is part of the MC_HRSR program.
  7              !
  8              ! Right-handed coordinates are assumed: X=down, Z=downstream, Y = (Z cross X)
  9              !
 10              ! Inputs are from common block in track_*.inc (except for fry):
 11              !  xs, ys, fry  are in cm.
 12              !  dxdzs, dydzs are unitless slopes (we say "radians" to contrast to "mr").
 13              ! Matrix Elements want meters and radians, so have to convert cm-->m.
 14              ! Output:
 15              !  delta_p is deviation from central momentum(unitless). Convert to % for output.
 16              !  delta_t, delta_phi are in "radians".
 17              !  y_tgt is in meters, convert to cm for output.
 18              !
 19              ! Author: D. Potterveld, ANL, 18-Mar-1993
 20              !
 21              ! Modification History:
 22 gaskelld 1.1 !
 23              ! 2-August-1995 (RMM) Hardcoded in an extra index on the coeff, expon, and
 24              !		      n_terms variables to specify HRSR. (HRSR = 3)
 25              !
 26              !  19-AUG-1993	(DHP) Modified to use COSY INFINITY reconstruction coefficients.
 27              C-______________________________________________________________________________
 28                   
 29              	implicit none
 30              
 31              	include '../track.inc'
 32              
 33              	integer*4 specnum
 34              	parameter (specnum = 3)			!this is the HRSR routine
 35              
 36              C Argument definitions.
 37              
 38              	real*8	delta_p,delta_t,delta_phi,y_tgt
 39              	real*8  fry			!vertical position at target (+y=down)
 40              
 41              C Cosy reconstruction matrix elements.
 42              
 43 gaskelld 1.1 	integer*4	max_elements
 44              	parameter	(max_elements = 1000)
 45              
 46              	integer*4	nspectr
 47              	parameter	(nspectr=4)
 48              
 49              	real*8		coeff(nspectr,4,max_elements)
 50              	integer*2	expon(nspectr,4,max_elements)
 51              C	integer*2	expon(nspectr,5,max_elements)
 52              	integer*4	n_terms(nspectr),max_order
 53              C	real*8		sum(4),hut(5),term
 54              	real*8		sum(4),hut(4),term
 55              
 56              C Misc. variables.
 57              
 58              	integer*4	i,j
 59              	integer*4	chan
 60              
 61              	logical		firsttime	/.true./
 62              	character*132	line
 63              
 64 gaskelld 1.1 C Functions.
 65              
 66              	logical locforunt
 67              
 68              C No amnesia, please...
 69              
 70              	save
 71              
 72              C ============================= Executable Code ================================
 73              
 74              C First time through, read in coefficients from data file.
 75              
 76              	if (firsttime) then
 77              	   if (.not.locforunt(chan)) stop 'MC_HRSR_RECON: No I/O channels!'
 78              	   open (unit=chan,status='old',readonly,file='hrsr/hrs_recon_cosy.dat')
 79              
 80              ! Skip past header.
 81              
 82              	   line = '!'
 83              	   do while (line(1:1).eq.'!')
 84              	      read (chan,1001) line
 85 gaskelld 1.1 	   enddo
 86              
 87              ! Read in coefficients and exponents.
 88              	   n_terms(specnum) = 0
 89              	   max_order = 0
 90              	   do while (line(1:4).ne.' ---')
 91              	      n_terms(specnum) = n_terms(specnum) + 1
 92              	      if (n_terms(specnum).gt.max_elements)
 93                   >	      stop 'WCRECON: too many COSY terms!'
 94              	      read (line,1200) (coeff(specnum,i,n_terms(specnum)),i=1,4),
 95                   >			       (expon(specnum,j,n_terms(specnum)),j=1,4)
 96              	      read (chan,1001) line
 97              	      max_order = max(max_order, expon(specnum,1,n_terms(specnum)) +
 98                   >				expon(specnum,2,n_terms(specnum)) +
 99                   >				expon(specnum,3,n_terms(specnum)) +
100                   >				expon(specnum,4,n_terms(specnum)))
101              C     >				expon(specnum,5,n_terms(specnum)))
102              	   enddo
103              	   write(6,*) 'HRSR:N_TERMS, MAX_ORDER = ',n_terms(specnum),max_order
104              	   close (unit=chan)
105              	   firsttime = .false.
106 gaskelld 1.1 	endif
107              
108              C Reset COSY sums.
109              
110              	do i = 1,4
111              	   sum(i) = 0.
112              	enddo
113              
114              C Convert hut quantities to right-handed coordinates, in meters and "radians".
115              	hut(1) = xs/100.		!cm --> m
116              	hut(2) = dxdzs			!slope ("radians")
117              	hut(3) = ys/100.		!cm --> m
118              	hut(4) = dydzs			!slope ("radians")
119              
120              C Compute COSY sums.
121              
122              	do i = 1,n_terms(specnum)
123              	   term = hut(1)**expon(specnum,1,i)*hut(2)**expon(specnum,2,i)
124                   >	        * hut(3)**expon(specnum,3,i)*hut(4)**expon(specnum,4,i)
125              	   sum(1) = sum(1) + term*coeff(specnum,1,i)
126              	   sum(2) = sum(2) + term*coeff(specnum,2,i)
127 gaskelld 1.1 	   sum(3) = sum(3) + term*coeff(specnum,3,i)
128              	   sum(4) = sum(4) + term*coeff(specnum,4,i)
129              	enddo
130                   
131              C Load output values.
132              
133              	delta_phi = sum(1)		!slope ("radians")
134              	y_tgt	  = sum(2)*100.		!m --> cm
135              	delta_t   = sum(3)		!slope ("radians")
136              	delta_p   = sum(4)*100.		!percent deviation
137              
138              	return
139              
140              C ============================ Format Statements ===============================
141              
142              1001	format(a)
143              1200	format(1x,4g16.9,1x,4i1)
144              
145                    END