**********************************
*           hcfit        
***************************
* reads in the beamposition
* from coda and harp runs
* and tries to calculate the 
* calibration factors
**********************************

v/del *

HARPFILE=/w/work506/e94104/espace/harp/harp_feb16_b.dat
CODAFILE=/w/work506/e94104/espace/harp/bpm_feb16_new.dat
psfile=/w/work506/e94104/espace/harp/hcfit_feb16_new.ps
sh rm [psfile]
for/file 99 [psfile]
meta 99 -111


*read HARPFILE 'enter harp-filename'
*read CODAFILE 'enter coda-filename'


v/read hrun,hcrun,x0,y0,hxar,hyar,hxbr,hybr,hxzr,hyzr,w1,w2 [HARPFILE]
v/read crun,bxa,dbxa,bya,dbya,bxb,dbxb,byb,dbyb,bxz,dbxz,byz,dbyz [CODAFILE]

*trace off 

*opt stat
opt grid
opt nfit
vscale hxar 1e-6 hxa
*vscale dhxar 1e-6 dhxa
vscale hyar 1e-6 hya
*vscale dhyar 1e-3 dhya
vscale hxbr 1e-6 hxb
*vscale dhxbr 1e-3 dhxb
vscale hybr 1e-6 hyb
*vscale dhybr 1e-3 dhyb
vscale hxzr 1e-6 hxz
vscale hyzr 1e-6 hyz
vscale hxzr 1e-3 hxz2
vscale hyzr 1e-3 hyz2

sigma nv=nco(hrun)
n= nv(1)
message number of harp runs [nv] 
vec/write [n]
sigma nv=nco(crun)
n=nv(1)
message number of coda runs [nv]
vec/write [n]

zone 2 2
set mtype 20
title ' ' 

vec/plot hrun%crun
vec/plot y0%x0
vec/plot hxz2%x0
vec/plot hyz2%y0

wait

*vec/write crun,hrun,distax,distay,distbx,distby,distzx,distzy ' ' (2f6.1,6(1x,*g12.5))
*message average difference in mm
*message bpmax , bpmay, bpmbx , bpmby, tar_x , tar_y   
*vec/write avdistax,avdistay,avdistbx,avdistby,avdistzx,avdistzy ' ' (6(1x,g12.*5))

zone 2 3
b0=1.0
v/cr par(2) R 0.0 [b0]
graph [n] bxa hxa 'W*'
v/fit bxa hxa dbxa p1 'QN' 2 par
atitle 'Harp AX vs Coda AX   '
message bpmax offset and slope
vec/write par ' ' (2g12.5)
sigma xmin=vmin(bxa)
sigma xmax=vmax(bxa)
m=xmin(1)
mm=xmax(1)
ym=par(1)+par(2)*[m]
ymm=par(1)+par(2)*[mm]
message [m] [mm] [ym] [ymm] 
dline [m] [mm] [ym] [ymm] 
null 0. 10. 0. 10. 'as'
t = a = $eval(par(1)) b = $eval(par(2)) 
text .5 10.5 [t] .2 0 l
wait
v/del par


v/cr par(2) R 0.0 [b0]
graph [n] bya hya 'W*'
v/fit bya hya dbya p1 'QN' 2 par
atitle 'Harp AY vs Coda AY   '
message bpmay offset and slope
vec/write par ' ' (2g12.5)
sigma xmin=vmin(bya)
sigma xmax=vmax(bya)
m=xmin(1)
mm=xmax(1)
ym=par(1)+par(2)*[m]
ymm=par(1)+par(2)*[mm]
*message [m] [mm] [ym] [ymm] 
dline [m] [mm] [ym] [ymm] 
null 0. 10. 0. 10. 'as'
t = a = $eval(par(1)) b = $eval(par(2)) 
text .5 10.5 [t] .2 0 l
wait
v/del par


v/cr par(2) R 0.0 [b0]
graph [n] bxb hxb 'W*'
v/fit bxb hxb dbxb p1 'QN' 2 par
atitle 'Harp BX vs Coda BX   '
message bpmbx offset and slope
vec/write par ' ' (2g12.5)
sigma xmin=vmin(bxb)
sigma xmax=vmax(bxb)
m=xmin(1)
mm=xmax(1)
ym=par(1)+par(2)*[m]
ymm=par(1)+par(2)*[mm]
*message [m] [mm] [ym] [ymm] 
dline [m] [mm] [ym] [ymm] 
null 0. 10. 0. 10. 'as'
t = a = $eval(par(1)) b = $eval(par(2)) 
text .5 10.5 [t] .2 0 l
wait
v/del par


v/cr par(2) R 0.0 [b0]
graph [n] byb hyb 'W*'
v/fit byb hyb dbyb p1 'QN' 2 par
atitle 'Harp BY vs Coda BY   '
message bpmby offset and slope
vec/write par ' ' (2g12.5)
sigma xmin=vmin(byb)
sigma xmax=vmax(byb)
m=xmin(1)
mm=xmax(1)
ym=par(1)+par(2)*[m]
ymm=par(1)+par(2)*[mm]
*message [m] [mm] [ym] [ymm] 
dline [m] [mm] [ym] [ymm] 
null 0. 10. 0. 10. 'as'
t = a = $eval(par(1)) b = $eval(par(2)) 
text .5 10.5 [t] .2 0 l
wait
v/del par


v/cr par(2) R 0.0 [b0]
graph [n] bxz hxz 'W*'
v/fit bxz hxz dbxz p1 'QN' 2 par
atitle 'Harp TarX vs Coda TarX   '
message bpmtarx offset and slope
vec/write par ' ' (2g12.5)
sigma xmin=vmin(bxz)
sigma xmax=vmax(bxz)
m=xmin(1)
mm=xmax(1)
ym=par(1)+par(2)*[m]
ymm=par(1)+par(2)*[mm]
*message [m] [mm] [ym] [ymm] 
dline [m] [mm] [ym] [ymm] 
null 0. 10. 0. 10. 'as'
t = a = $eval(par(1)) b = $eval(par(2)) 
text .5 10.5 [t] .2 0 l
wait
v/del par



v/cr par(2) R 0.0 [b0]
graph [n] byz hyz 'W*'
v/fit byz hyz dbyz p1 'QN' 2 par
atitle 'Harp TarY vs Coda TarY   '
message bpmtary offset and slope
vec/write par ' ' (2g12.5)
sigma xmin=vmin(byz)
sigma xmax=vmax(byz)
m=xmin(1)
mm=xmax(1)
ym=par(1)+par(2)*[m]
ymm=par(1)+par(2)*[mm]
*message [m] [mm] [ym] [ymm] 
dline [m] [mm] [ym] [ymm] 
null 0. 10. 0. 10. 'as'
t = a = $eval(par(1)) b = $eval(par(2)) 
text .5 10.5 [t] .2 0 l
wait
v/del par

close 99