;+
; NAME:
;   plot_sxt_vig
; PURPOSE:
;   Make the plots for SXT Cal Note 37 (Vignette function)
; CALLING SEQUENCE:
;   plot_sxt_vig,/high
;   plot_sxt_vig,/high,/ps
;   plot_sxt_vig,/low
;   plot_sxt_vig,/low,/ps
;   plot_sxt_vig,/low,/ps,path=path  ; If data files not on current directory
;   plot_sxt_vig,/low,/ps,/noplot    ; Create plot file but don't print
; INPUT KEYWORDS:
;   high	- If set, compute the high-energy case
;   low		- If set, compute the low-energy case
;   ps		- If set, plot the results.  Also create two files:
;			cal_note37_hi.ps and cal_note37_lo.ps
;   path	- To specify the directory containing the data files.
;
; PROCEDURE:
;   This is a very specialized routine to generate the final plots
;   of the Vignette study that was conducted by R. A. Fuller.  His
;   subroutines are included in this file.
;   
;   These routine require several data files to work.  They are available,
;   in the on-line directories, but the user will have to copy the compressed 
;   Tar file to his or her own scratch area.  The Tar file must be uncompressed 
;   and the files must be extracted.
;
; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;   Instructions to restore the data file:
;
;   The data files are contained in a compressed Tar file called
;   fuller_vignette.tar.Z on the $DIR_SXT_SENSITIVE directory.
;
;   To run this program you must:
;   1. Copy the fuller_vignette.tar.Z file from $DIR_SXT_SENSITIVE
;      to a local directory. It might be best to create a new sub-directory.
;      Most of files will be names 'wc_*.dat' or 'wc_*.pc' and there is
;      a README file. On a Unix machine use the following commands:
;      % mkdir junk		# Optional -- create a sub-directory called junk
;      % cd junk		# Move to the sub-directory
;      % cp $DIR_SXT_SENSITIVE/fuller_vignette.tar.Z .
;   2. Make sure you are on the sub-directory which contains the compressed tar file.
;      % ls -l fuller_vignette*
;   3. At the Unix prompt type:
;      % uncompress fuller_vignette.tar.Z
;   4. To unpack the Tar file, at the Unix prompt type:
;      % tar xvf fuller_vignette.tar
;   5. Run this program in IDL.  If you are not running from the same
;      directory as the data files, use the path=path to specify the
;      location of the files.
;      % idl
;      IDL> plot_sxt_vig,/lo	; For low-energy plot (see calling sequence)
; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;
;
;   In these programs the centroid positions and the slopes have been
;   fixed to the best values.  The routine will fit the data point to
;   these curves and compute the Chi square values.
;
;   The /special switch eliminates 3 data points that have partial annuli
;   yet are within the nominal "inner-region".  This was a Fuller fix for
;   the fact that the inner- and outer-region cones are not concentric.
;   That is, the /special switch removes points that are really in the
;   outer region.
;
; RESTRICTIONS/SIDE EFFECTS:
;   Will try to set the plot display to 'X' for diagnostic plots.
;   Not guaranteed on VMS.
; MODIFICATION HISTORY:
;   2-Feb-94, J. R. Lemen, Written
;-
;
;  *** plot_sxt_vig appears as the last routine ***
;
;

function average, x

	return, total(x)/N_ELEMENTS(x) 
end;

pro eff_area_lo, ps=ps, XC = XC, YC = YC, Rcut = Rcut, nsig=nsig, $
	special=special, oXCo=xXCo, oYCo=oYCo, Rmax=Rmax, slope=slope, $
	nofix=nofix, path=path, noplot=noplot, fileout=fileout, qtest=qtest
;+
; MODIFICATION HISTORY:
; 21-jul-93, R. A. Fuller, Written
; 29-jan-94, RAF + JRL,  Fixed defaults for least scatter
;			 run with /special to get best case
; 31-jan-94, JRL, Renamed to eff_area_lo; Pass in the path for the data files.
;  2-feb-94, JRL, Added noplot and fileout keywords
;-


   if n_elements(Rmax) eq 0 then Rmax = 24.		; 29-jan-94
   if n_elements(Rcut) eq 0 then Rcut = 17.1		; 29-jan-94

   if N_ELEMENTS(Rcut) then if (Rcut LT 30.0) then $
			Rcut = Rcut/gt_pix_size()*60.
   if N_ELEMENTS(Rmax) then if (Rmax LT 100.0) then $
			Rmax = Rmax/gt_pix_size()*60.
   if NOT(N_ELEMENTS(XC)) then XC = 515.;
   if NOT(N_ELEMENTS(YC)) then YC = 633.;
   if NOT(N_ELEMENTS(Rcut)) then Rcut = 18.5*60/gt_pix_size();
;   if NOT(N_ELEMENTS(nsig)) then nsig = 1.0;
   if NOT(N_ELEMENTS(nsig)) then nsig = 3.0;		; 29-jan-94
   if NOT(N_ELEMENTS(oXCo)) then oXCo = 533.;
   if NOT(N_ELEMENTS(oYCo)) then oYCo = 599.;



   if n_elements(path) eq 0 then begin			; 31-jan-94
	path = '/3p/fuller/data'
        path2 = concat_dir( path, 'pc')
   endif else path2 = path

   fn1 = concat_dir(path,'ws_al_7x7y_out_01rs.dat')
   fn2 = concat_dir(path,'ws_al_7x7y_in_01rs.dat')
   fn3 = concat_dir(path,'ws_c_7x7y_best_01rs.dat')
   fn4 = concat_dir(path,'ws_al_5x5y_in_01.dat')
   fn5 = concat_dir(path,'ws_al_4x4y_out_01.dat')
  

;========================================================================


   print, format = '(/,/," Reading file :", A, /)', fn1
   restore, file=fn1, /verbose
   break_file, fn1, disk, dir, buf, ext		; 2-feb-94

;   len = strpos(fn1, '.dat') - strpos(fn1, 'data')-5
;   buf = strmid(fn1, strpos(fn1, 'data')+5, len)
   if (strpos(buf, 'rs') GT 0) then buf = strmid(buf, 0, strpos(buf, 'rs'))

   d_al_7x7y_out = gdataorig.bssum;   
   x_al_7x7y_out = gdataorig.xcenter;
   y_al_7x7y_out = gdataorig.ycenter;


   corr4flux, buf, x_al_7x7y_out, y_al_7x7y_out, d_al_7x7y_out, nsig=nsig, path=path2

   if N_ELEMENTS(special) then begin
	ii = indgen(46);
        ii(45) = 48;
	ii(44) = 47;
        d_al_7x7y_out = d_al_7x7y_out(ii);
        x_al_7x7y_out = x_al_7x7y_out(ii);
        y_al_7x7y_out = y_al_7x7y_out(ii);
   endif

   errorBars = 1./sqrt(d_al_7x7y_out);
  
   b_al_7x7y_out = 8227.25;
   d_al_7x7y_out = d_al_7x7y_out/b_al_7x7y_out;
   n_al_7x7y_out = N_ELEMENTS(d_al_7x7y_out);

   data = d_al_7x7y_out;
   x = x_al_7x7y_out;
   y = y_al_7x7y_out;


;========================================================================



   print, format = '(/,/," Reading file :", A, /)', fn2
   restore, file=fn2, /verbose
   break_file, fn2, disk, dir, buf, ext		; 2-feb-94

;   len = strpos(fn2, '.dat') - strpos(fn2, 'data')-5
;   buf = strmid(fn2, strpos(fn2, 'data')+5, len)
   if (strpos(buf, 'rs') GT 0) then buf = strmid(buf, 0, strpos(buf, 'rs'))

   d_al_7x7y_in = gdataorig.bssum;   
   x_al_7x7y_in = gdataorig.xcenter;
   y_al_7x7y_in = gdataorig.ycenter;


   corr4flux, buf, x_al_7x7y_in, y_al_7x7y_in, d_al_7x7y_in, nsig=nsig, path=path2

   if N_ELEMENTS(special) then begin
	ii = indgen(46);
        ii(45) = 48;
	ii(44) = 47;
        d_al_7x7y_in = d_al_7x7y_in(ii);
        x_al_7x7y_in = x_al_7x7y_in(ii);
        y_al_7x7y_in = y_al_7x7y_in(ii);
   endif

   errorBars = [errorBars,1./sqrt(d_al_7x7y_in)];

   b_al_7x7y_in = 7591.33;
   d_al_7x7y_in = d_al_7x7y_in/b_al_7x7y_in;
   n_al_7x7y_in = N_ELEMENTS(d_al_7x7y_in);

   data = [data,d_al_7x7y_in];
   x = [x,x_al_7x7y_in];
   y = [y,y_al_7x7y_in];


;========================================================================



   print, format = '(/,/," Reading file :", A, /)', fn3
   restore, file=fn3, /verbose
   break_file, fn3, disk, dir, buf, ext			; 2-feb-94

;   len = strpos(fn3, '.dat') - strpos(fn3, 'data')-5
;   buf = strmid(fn3, strpos(fn3, 'data')+5, len)
   if (strpos(buf, 'rs') GT 0) then buf = strmid(buf, 0, strpos(buf, 'rs'))

   d_c_7x7y_best = gdataorig.bssum;   
   x_c_7x7y_best = gdataorig.xcenter;
   y_c_7x7y_best = gdataorig.ycenter;


   corr4flux, buf, x_c_7x7y_best, y_c_7x7y_best, d_c_7x7y_best, nsig=nsig, path=path2

   if N_ELEMENTS(special) then begin
	ii = indgen(46);
        ii(45) = 48;
  	 ii(44) = 47;
        d_c_7x7y_best = d_c_7x7y_best(ii);
        x_c_7x7y_best = x_c_7x7y_best(ii);
        y_c_7x7y_best = y_c_7x7y_best(ii);
  endif

   errorBars = [errorBars,1./sqrt(d_c_7x7y_best)];

   b_c_7x7y_best = 5964.49;
   d_c_7x7y_best = d_c_7x7y_best/b_c_7x7y_best;
   n_c_7x7y_best = N_ELEMENTS(d_c_7x7y_best);

   data = [data,d_c_7x7y_best];
   x = [x,x_c_7x7y_best];
   y = [y,y_c_7x7y_best];



;========================================================================

   print, format = '(/,/," Reading file :", A, /)', fn4
   restore, file=fn4, /verbose
   break_file, fn4, disk, dir, buf, ext			; 2-feb-94

;   len = strpos(fn4, '.dat') - strpos(fn4, 'data')-5
;   buf = strmid(fn4, strpos(fn4, 'data')+5, len)
   if (strpos(buf, 'rs') GT 0) then buf = strmid(buf, 0, strpos(buf, 'rs'))

   d_al_5x5y_in = gdata.bssum;   
   x_al_5x5y_in = gdata.xcenter;
   y_al_5x5y_in = gdata.ycenter;


   corr4flux, buf, x_al_5x5y_in, y_al_5x5y_in, d_al_5x5y_in, nsig=nsig, path=path2

   errorBars = [errorBars,1./sqrt(d_al_5x5y_in)];

   b_al_5x5y_in = 7683.9;
   d_al_5x5y_in = d_al_5x5y_in/b_al_5x5y_in;
   n_al_5x5y_in = N_ELEMENTS(d_al_5x5y_in);

   data = [data,d_al_5x5y_in];
   x = [x,x_al_5x5y_in];
   y = [y,y_al_5x5y_in];


;========================================================================



   print, format = '(/,/," Reading file :", A, /)', fn5
   restore, file=fn5, /verbose
   break_file, fn5, disk, dir, buf, ext 			; 2-feb-94

;   len = strpos(fn5, '.dat') - strpos(fn5, 'data')-5
;   buf = strmid(fn5, strpos(fn5, 'data')+5, len)
   if (strpos(buf, 'rs') GT 0) then buf = strmid(buf, 0, strpos(buf, 'rs'))

   d_al_4x4y_out = gdata.bssum;   
   x_al_4x4y_out = gdata.xcenter;
   y_al_4x4y_out = gdata.ycenter;


   corr4flux, buf, x_al_4x4y_out, y_al_4x4y_out, d_al_4x4y_out, nsig=nsig, path=path2

   errorBars = [errorBars,1./sqrt(d_al_4x4y_out)];

   b_al_4x4y_out = 8046.8;
   d_al_4x4y_out = d_al_4x4y_out/b_al_4x4y_out;
   n_al_4x4y_out = N_ELEMENTS(d_al_4x4y_out);

   data = [data,d_al_4x4y_out];
   x = [x,x_al_4x4y_out];
   y = [y,y_al_4x4y_out];


;========================================================================
    genb = (b_al_7x7y_out+b_al_7x7y_in+b_c_7x7y_best+b_al_5x5y_in + $
		b_al_4x4y_out)/5.;
;    genb = (b_al_7x7y_out+b_al_7x7y_in+b_al_5x5y_in + $
;		b_al_4x4y_out)/4.;
;    genb = (b_al_7x7y_out+b_al_7x7y_in+b_c_7x7y_best)/3.;

;========================================================================


   R = sqrt((x-XC)^2 + (y-YC)^2)	; Get R
   sortArr = sort(R);

   R = R(sortArr);
   ii = where(R LT Rcut);
   R = R(ii);
   sdatin = data(sortArr);
   sdatin = sdatin(ii);
   eBin = errorBars(sortArr);
   eBin = eBin(ii);

   pp = fltarr(2);

;   if NOT(N_ELEMENTS(slope)) then slope = -0.000490;
;   if NOT(N_ELEMENTS(slope)) then slope = -0.0005085	; 29-jan-94 (0.01244 per arc-min)
   if NOT(N_ELEMENTS(slope)) then slope = -0.000510829	; 29-jan-94 (0.01244 per arc-min)

   pp(0) = 1.0;
   pp(1) = slope*pp(0);
;
;  1/sd^2 = Instrumental Weighting
;  1/y = Statistical Weighting
;  1.0 = No weighting
;
   WF = replicate(1., N_ELEMENTS(sdatin));
   ii = where(sdatin gt 0.) & WF(ii) = 1./sdatin(ii)

   dFit = curvefit(R(0:*), sdatin(0:*), WF(0:*),  pp, sig, funct="line")
   if NOT(N_ELEMENTS(nofix)) then begin
	;
	;  Revert to desired slope
	;
	   pp(1) = slope;
	   dFit = pp(0) + pp(1)*R;
   endif;
;
;
   redchi = genb*total( WF*(sdatin-dFit)^2)/(N_ELEMENTS(sdatin) -2);

   corr = abs(correlate(R,sdatin));


   xout = x;
   yout = y;

   Rout = sqrt((xout-oXCo)^2 + (yout-oYCo)^2)	; Get Rout
   sortArr = sort(Rout);
   Rout = Rout(sortArr);
   sdatout = data(sortArr);
   eBout = errorBars(sortArr);
   ii = where(Rout GE Rcut);
   Rout = Rout(ii);
   sdatout = sdatout(ii);
   eBout = eBout(ii);
   jj = where(Rout LE Rmax);
   Rout = Rout(jj);
   sdatout = sdatout(jj);
   eBout = eBout(jj);

   ppout = fltarr(2);

   ppout(0) = max(sdatout);
;   ppout(1) = -0.00180894*pp(0);
;   ppout(1) = -0.0017472*pp(0);	29-jan-94 (0.04274 /arc-min)
   ppout(1) = -0.00175520*pp(0);	29-jan-94 (0.04274 /arc-min)

;
;  1/sd^2 = Instrumental Weighting
;  1/y = Statistical Weighting
;  1.0 = No weighting
;
   WF = replicate(1., N_ELEMENTS(sdatout));
   ii = where(sdatout gt 0.) & WF(ii) = 1./sdatout(ii)

   dFitout = curvefit(Rout(0:*), sdatout(0:*), WF(0:*),  ppout, sigout, funct="line")

   redchiout = genb*total( WF*(sdatout-dFitout)^2)/(N_ELEMENTS(sdatout) -2);

   corrout = abs(correlate(Rout,sdatout));


   print,'Inner-offset = ',pp(0)
   b = pp(0);
   pp(0) = 1.0;
   pp(1) = (pp(1)/b)*(60./gt_pix_size())
   ppout(0) = ppout(0)/b;
   ppout(1) = (ppout(1)/b)*(60./gt_pix_size())
   Rcrit = (ppout(0) - pp(0))/(pp(1) - ppout(1));

   Ri = [0, Rcrit];
   dati = pp(0) + Ri*pp(1);
   Ro = [Rcrit , 25];
   dato = ppout(0) + Ro*ppout(1);

   R = R*gt_pix_size()/60;
   Rout = Rout*gt_pix_size()/60;
   RR = [R,Rout];
   sdattot = [sdatin,sdatout];
   eBtot = [eBin,eBout];
   WF = replicate(1., N_ELEMENTS(sdattot));
   ii = where(sdattot gt 0.) & WF(ii) = 1./sdattot(ii)
   redchitot=genb*total(WF*(sdattot-[dFit,dFitout])^2)/(N_ELEMENTS(sdattot)-2);

   sdattot = sdattot/b;

   if N_ELEMENTS(ps) then begin
	print, '--------------- Making PostScript Plot ------------------'
	set_plot,'ps'
	if n_elements(fileout) eq 0 then fileout = 'idl.ps'
	device, filename=fileout
	device, xsize=25.4, ysize=19.05
	device, /landscape
	device, /HELVETICA, /BOLD
	device, font_size = 14
   endif else begin
     set_plot,'X' & wdef,0,512
   endelse

	ang = findgen(16)*(!pi * 2/16.0);
   	usersym, 0.4*cos(ang), 0.4*sin(ang), /fill
   	plot, RR, sdattot, psym=8, xrange = [0,30], $
      		title =  'Effective Area of SXT vs Off Axis Angle', $
      		xtitle = 'Off Axis Angle ('')', $
		ytitle = 'Relative Effective Area',$
      		/data, yrange = [0,1],charthick = 2.0
   	errplot, RR, sdattot+eBtot, sdattot-eBtot
        oplot, Ri,dati
   	oplot, Ro,dato, linestyle = 2
        oplot, [0,40], 1-0.012*[0,40],linestyle=3
        oplot, [0,40], 1.648-0.0428*[0,40],linestyle=4

   	if (N_ELEMENTS(msg)) then xyouts, 0.02, 0.0, msg, /normal, charsize = 0.5
   
   	buffer = string(format='("Maximum Effective Area (Center) @:", 2F8.2)', $
		XC, YC);
   	xyouts, 0.15, 0.50, buffer, /normal, charsize = 0.75
   
   	buffer = string(format='("Maximum Effective Area (Outer)@:", 2F8.2)', $
		oXCo, oYCo);
   	xyouts, 0.15, 0.47, buffer, /normal, charsize = 0.75

   	buffer = string($
               format='("Inner-Fit to line : Yi=",F8.5," ",F8.5,"*Ri")', pp(0),pp(1));
   	xyouts, 0.15, 0.40, buffer, /normal, charsize = 0.75

   	buffer = string($

                format='("Reduced ChiSquared = ",F8.5)', redchi);
   	xyouts, 0.15, 0.37, buffer, /normal, charsize = 0.75

   	buffer = string($

                format='("Correlation Coefficient = ",F8.5)', corr);
   	xyouts, 0.15, 0.34, buffer, /normal, charsize = 0.75


   	buffer = string(format='("---.---.---.   0.012 slope line")');
   	xyouts, 0.65, 0.37, buffer, /normal, charsize = 0.75

        buffer = string(format='("--...--...--...   0.0428 slope line.")');
   	xyouts, 0.65, 0.34, buffer, /normal, charsize = 0.75


   	buffer = string($
                format='("Correction Tolerance = ",F4.1," sigma")', nsig);
   	xyouts, 0.65, 0.28, buffer, /normal, charsize = 0.75


        if N_ELEMENTS(Rmax) then begin
	   	buffer = string($
	                format='("Maxium Radius = ",F5.2,"''")', $
			Rmax*gt_pix_size()/60.);
   		xyouts, 0.65, 0.15, buffer, /normal, charsize = 0.75
        endif;


        buffer = string($
               format='("Outer-Fit to line : Yo=",F8.5," ",F8.5,"*Ro")', ppout(0),ppout(1));
   	xyouts, 0.15, 0.23, buffer, /normal, charsize = 0.75

   	buffer = string($

                format='("Reduced ChiSquared = ",F8.5)', redchiout);
   	xyouts, 0.15, 0.20, buffer, /normal, charsize = 0.75

   	buffer = string($

                format='("Correlation Coefficient = ",F8.5)', corrout);
   	xyouts, 0.15, 0.17, buffer, /normal, charsize = 0.75


   	buffer = string($

                format='("Combined Reduced ChiSquared = ",F8.5)', redchitot);
   	xyouts, 0.15, 0.12, buffer, /normal, charsize = 0.75


   	buffer = string($

                format='("Critical Radius = ",F5.2, "''")', Rcrit);
   	xyouts, Rcrit+0.05, dati(1)+0.05, buffer, /data, charsize = 0.75

	buffer = string($
                format='("Cutoff Radius = ",F5.2, "''")', Rcut*gt_pix_size()/60.);
   	xyouts, Rcut*gt_pix_size()/60.+0.05, dati(1)+0.08, buffer, /data, charsize = 0.75

        xyouts, 0.02, 0.0, 'al_7x7y_in, al_7x7y_out, c01_7x7y_best' + $
		', al_5x5y_in, al_4x4y_out', /normal, charsize = 0.5
;        xyouts, 0.02, 0.0, 'al_7x7y_in, al_7x7y_out' + $
;		', al_5x5y_in, al_4x4y_out', /normal, charsize = 0.5

   	plottime

   if N_ELEMENTS(ps) then begin
	device, font_size = 12

	fff = fstat(!d.unit)
	message,'Plot file written: '+fff.name,/info
	device, /close
	if not keyword_set(noplot) then pprint, fileout		;, /delete

	set_plot,'x'
   endif


   print, N_ELEMENTS(data), N_ELEMENTS(sdatin),N_ELEMENTS(sdatout);
   print, N_ELEMENTS(sdatin) + N_ELEMENTS(sdatout);

   if keyword_set(qtest) then begin
      ans='' & read,'Hit CR to continue',ans
   endif


   return;
end;

pro eff_area_hi, ps=ps, XC = XC, YC = YC, oXCo=oXCo, oYCo=oYCo, $
          Rcut = Rcut, nsig=nsig, special=special, slope=slope, $
 	  reshape=reshape, Rmax = Rmax, nofix = nofix, 		$
	  path=path, noplot=noplot, fileout=fileout, qtest=qtest
;+
; MODIFICATION HISTORY
; 19-aug-93, R. A. Fuller, Written
; 29-jan-94, RAF+JRL, Modified defaults to give best results (slope and Rmax)
; 31-jan-94, JRL, Renamed to eff_area_hi; Pass in the path for the data files
;  2-feb-94, JRL, Added noplot and fileout keywords
;-

   if N_elements(Rmax) eq 0 then Rmax = 24.0		; 29-jan-94

   if N_ELEMENTS(Rcut) then if (Rcut LT 30.0) then $
			Rcut = Rcut/gt_pix_size()*60.
   if N_ELEMENTS(Rmax) then if (Rmax LT 100.0) then $
			Rmax = Rmax/gt_pix_size()*60.
   if NOT(N_ELEMENTS(XC)) then XC = 515.;
   if NOT(N_ELEMENTS(YC)) then YC = 629.;
   if NOT(N_ELEMENTS(Rcut)) then Rcut = 18.4*60./gt_pix_size();
   if NOT(N_ELEMENTS(nsig)) then nsig = 1.0;
   if NOT(N_ELEMENTS(oXCo)) then oXCo = 533.;
   if NOT(N_ELEMENTS(oYCo)) then oYCo = 599.;		; 2-feb-94 (Changed from 598)


;   fn = pickfile(path = '/3p/fuller/data', filter = '*.dat', $
;              /NOCONFIRM, /READ);
   if n_elements(path) eq 0 then begin				; 31-jan-94
	path = '/3p/fuller/data'
   	path2 = concat_dir( path, 'pc')
   endif else path2 = path
   fn = concat_dir(path,  'ws_ag_7x7y_best_01rs.dat')		; 31-jan-94
   print, format = '(/,/," Reading file :", A, /)', fn
   if (strlen(fn) GT 19) then restore, file=fn, /verbose $
   else return

   if N_ELEMENTS(gdataorig) then begin
	   sdat = gdataorig.bssum;   
	   x = gdataorig.xcenter;
	   y = gdataorig.ycenter;
   endif else begin
	   sdat = gdata.bssum;   
	   x = gdata.xcenter;
	   y = gdata.ycenter;
   endelse

   break_file, fn, disk, dir, buf, ext 				; 2-feb-94

;   len = strpos(fn, '.dat') - strpos(fn, 'data')-5
;   buf = strmid(fn, strpos(fn, 'data')+5, len)
   if (strpos(buf, 'rs') GT 0) then buf = strmid(buf, 0, strpos(buf, 'rs'))

   corr4flux, buf, x, y, sdat, nsig=nsig, path = path2		; 31-jan-94

   if N_ELEMENTS(special) then begin
	ii = indgen(46);
        ii(45) = 48;
	ii(44) = 47;
        sdat = sdat(ii);
        x = x(ii);
        y = y(ii);
   endif

   if N_ELEMENTS(reshape) then begin
        center_points,sdat
        center_points,x
        center_points,y
   endif

   if N_ELEMENTS(Rmax) then begin
	   R = sqrt((x-XC)^2 + (y-YC)^2)	; Get R
	   jj = where(R LE Rmax);
	   sdat = sdat(jj);
	   x = x(jj);
	   y = y(jj);
   endif;


   n = N_ELEMENTS(sdat);

   errorBars = [1./sqrt(sdat)]


;   if NOT(N_ELEMENTS(slope)) then slope = -0.0004840;
   if NOT(N_ELEMENTS(slope)) then slope = -0.0009545	; 29-jan-94
   pp = fltarr(2);
   pp(0) = max(sdat);
   pp(1) = slope*pp(0);

   R = sqrt((x-XC)^2 + (y-YC)^2)	; Get R
   sortArr = sort(R);
   R = R(sortArr);
   ii = where(R LT Rcut);
   if (ii(0) GT -1) then begin
	   R = R(ii);
	   sdatin = sdat(sortArr);
	   sdatin = sdatin(ii);
	   eBin = errorBars(sortArr);
	   eBin = eBin(ii);
   endif;


   if (ii(0) GT -1) then begin
	;
	;  1/sd^2 = Instrumental Weighting
	;  1/y = Statistical Weighting
	;  1.0 = No weighting
	;
	   WF = replicate(1., N_ELEMENTS(sdatin));
	   ii = where(sdatin gt 0.) & WF(ii) = 1./sdatin(ii)

	   dFit = curvefit(R(0:*), sdatin(0:*), WF(0:*),pp, sig, funct="line")
           if NOT(N_ELEMENTS(nofix)) then begin
		   ; Revert to previous slope but use the fit intercept
		   pp(1) = pp(0) * slope;
                   dFit = pp(0) + pp(1)*R;
           endif;

	   redchi = total( WF*(sdatin-dFit)^2)/(N_ELEMENTS(sdatin) -2);
	   Rbar = average(R);
	   ybar = average(sdatin);
	   sdr = total( (R-Rbar)^2)
	   sdy = total( (sdatin-ybar)^2)
	   sdRy = total( (R-Rbar)*(sdatin-ybar))

	   corr = abs(sdRy/sqrt(sdR*sdy))
   endif else begin
	redchi = 0.0;
        corr = 1.0;
   endelse;


   ppout = fltarr(2);
   ppout(0) = 1.54064*pp(0);
   ppout(1) = -0.00178114*pp(0);

   Rout = sqrt((x-oXCo)^2 + (y-oYCo)^2)	; Get Rout
   sortArr = sort(Rout);
   Rout = Rout(sortArr);
   jj = where(Rout GE Rcut);
   if (jj(0) GT -1) then begin
	   Rout = Rout(jj);
	   sdatout = sdat(sortArr);
	   sdatout = sdatout(jj);
	   eBout = errorBars(sortArr);
	   eBout = eBout(jj);
   endif;

   if (jj(0) GT -1) then begin
	;
	;  1/sd^2 = Instrumental Weighting
	;  1/y = Statistical Weighting
	;  1.0 = No weighting
	;
	   WF = replicate(1., N_ELEMENTS(sdatout));
	   ii = where(sdatout gt 0.) & WF(ii) = 1./sdatout(ii)

	   dFitout = curvefit(Rout(0:*), sdatout(0:*), $
				WF(0:*),  ppout, sigout, funct="line")
   
   	    redchiout=total( WF*(sdatout-dFitout)^2)/(N_ELEMENTS(sdatout) -2);

	   Rbar = average(Rout);
	   ybar = average(sdatout);
	   sdr = total( (Rout-Rbar)^2)
	   sdy = total( (sdatout-ybar)^2)
	   sdRy = total( (Rout-Rbar)*(sdatout-ybar))

	   corrout = abs(sdRy/sqrt(sdR*sdy))
   endif else begin
	redchiout = 0.0;
	corrout = 1.0;
   endelse;


   b = pp(0);
   pp(0) = 1.0;
   pp(1) = (pp(1)/b)*(60./gt_pix_size())
   ppout(0) = ppout(0)/b;
   ppout(1) = (ppout(1)/b)*(60./gt_pix_size())
   Rcrit = (ppout(0) - pp(0))/(pp(1) - ppout(1));

   Ri = [0, Rcrit];
   dati = pp(0) + Ri*pp(1);
   Ro = [Rcrit , 40];
   dato = ppout(0) + Ro*ppout(1);

   R = R*gt_pix_size()/60;
   Rout = Rout*gt_pix_size()/60;
   RR = [R,Rout];
   if (ii(0) GT -1) AND (jj(0) GT -1) then begin
	   sdattot = [sdatin,sdatout];
	   eBtot = [eBin,eBout];
           dFtot = [dFit,dFitout];
   endif else if (ii(0) GT -1) then begin
	   sdattot = sdatin;
	   eBtot = eBin;
           dFtot = dFit;
   endif else begin
	   sdattot = sdatout;
	   eBtot = eBout;
           dFtot = dFitout;
   endelse;

   WF = replicate(1., N_ELEMENTS(sdattot));
   ii = where(sdattot gt 0.) & WF(ii) = 1./sdattot(ii)
   redchitot=total(WF*(sdattot-dFtot)^2)/(N_ELEMENTS(sdattot)-2);

   sdattot = sdattot/b;

   if N_ELEMENTS(ps) then begin
	print, '--------------- Making PostScript Plot ------------------'
	set_plot,'ps'
	if n_elements(fileout) eq 0 then fileout = 'idl.ps'
	device, filename= fileout
	device, xsize=25.4, ysize=19.05
	device, /landscape
	device, /HELVETICA, /BOLD
	device, font_size = 14
   endif else begin
        set_plot,'X' & wdef,0,512
   endelse

	ang = findgen(16)*(!pi * 2/16.0);
   	usersym, 0.4*cos(ang), 0.4*sin(ang), /fill
;        if N_ELEMENTS(Rmax) then xmax = Rmax*gt_pix_size()/60. $
;	else xmax=max(RR);
   	plot, RR, sdattot, psym=8, xrange = [0,30], $
      		title =  'Effective Area of SXT vs Off Axis Angle', $
      		xtitle = 'Off Axis Angle ('')', $
		ytitle = 'Relative Effective Area',$
      		/data, yrange = [0,1],charthick = 2.0
   	errplot, RR, sdattot+eBtot, sdattot-eBtot
        oplot, Ri,dati
   	oplot, Ro,dato, linestyle = 2
        oplot, [0,40], 1-0.012*[0,40],linestyle=3
        oplot, [0,40], 1.648-0.0428*[0,40],linestyle=4

   	if (N_ELEMENTS(msg)) then xyouts, 0.02, 0.0, msg, /normal, charsize = 0.5
   
   	buffer = string(format='("Maximum Effective Area (Center) @:", 2F8.2, " Flight")', $
		XC, YC);
   	xyouts, 0.15, 0.50, buffer, /normal, charsize = 0.75
   
   	buffer = string(format='("Maximum Effective Area (Outer)@:",2F8.2," Flight")', $
		oXCo, oYCo);
   	xyouts, 0.15, 0.47, buffer, /normal, charsize = 0.75

   	buffer = string($
               format='("Inner-Fit to line : Yi=",F8.5," ",F8.5,"*Ri")', pp(0),pp(1));
   	xyouts, 0.15, 0.43, buffer, /normal, charsize = 0.75

   	buffer = string($

                format='("Reduced ChiSquared = ",F8.5)', redchi);
   	xyouts, 0.15, 0.40, buffer, /normal, charsize = 0.75

   	buffer = string($

                format='("Correlation Coefficient = ",F8.5)', corr);
   	xyouts, 0.15, 0.37, buffer, /normal, charsize = 0.75


   	buffer = string(format='("---.---.---.   0.012 slope line")');
   	xyouts, 0.15, 0.30, buffer, /normal, charsize = 0.75

        buffer = string(format='("--...--...--...   0.0428 slope line.")');
   	xyouts, 0.15, 0.28, buffer, /normal, charsize = 0.75


   	buffer = string($
                format='("Correction Tolerance = ",F4.1," sigma")', nsig);
   	xyouts, 0.65, 0.30, buffer, /normal, charsize = 0.75


        if N_ELEMENTS(Rmax) then begin
	   	buffer = string($
	                format='("Maxium Radius = ",F5.2,"''")', $
			Rmax*gt_pix_size()/60.);
   		xyouts, 0.65, 0.20, buffer, /normal, charsize = 0.75
        endif;


        buffer = string($
               format='("Outer-Fit to line : Yo=",F8.5," ",F8.5,"*Ro")', ppout(0),ppout(1));
   	xyouts, 0.15, 0.23, buffer, /normal, charsize = 0.75

   	buffer = string($

                format='("Reduced ChiSquared = ",F8.5)', redchiout);
   	xyouts, 0.15, 0.20, buffer, /normal, charsize = 0.75

   	buffer = string($

                format='("Correlation Coefficient = ",F8.5)', corrout);
   	xyouts, 0.15, 0.17, buffer, /normal, charsize = 0.75


   	buffer = string($

                format='("Combined Reduced ChiSquared = ",F8.5)', redchitot);
   	xyouts, 0.15, 0.12, buffer, /normal, charsize = 0.75


   	buffer = string($

                format='("Critical Radius = ",F5.2, "''")', Rcrit);
   	xyouts, Rcrit+0.05, dati(1)+0.05, buffer, /data, charsize = 0.75

	buffer = string($

                format='("Cutoff Radius = ",F5.2, "''")', Rcut*gt_pix_size()/60.);
   	xyouts, Rcut*gt_pix_size()/60.+0.05, dati(1)+0.08, buffer, /data, charsize = 0.75

        xyouts, 0.02, 0.0, fn, /normal, charsize = 0.5

   	plottime

   if N_ELEMENTS(ps) then begin
	device, font_size = 12

	fff = fstat(!d.unit)
	message,'Plot file written: '+fff.name,/info
	device, /close
	if not keyword_set(noplot) then pprint, fileout		;, /delete

	set_plot,'x'
   endif


   print, N_ELEMENTS(sdat), N_ELEMENTS(sdatin),N_ELEMENTS(sdatout);
   print, N_ELEMENTS(sdatin) + N_ELEMENTS(sdatout);

   if keyword_set(qtest) then begin
      ans='' & read,'Hit CR to continue',ans
   endif

   return
end;

;	(12-aug-92)
;
;+
;  Name:
;	CORR4FLUX.PRO
;
;  Purpose:
;	Analyze SXT Calibration Data 
;
;
;
;
;  Input:
;       Interactively prompts the user for information
;
;  Optional Input:
;
;      BUF   This parameter is a string with the name of the logical
;            filename to be extracted and analyzed
;
;      SAVE  If this parameter is non-zero, the data is saved
;
;
;
;
;   Written :
;		RAF   Aug 12, 1992 
;	31-jan-94, J. R. Lemen, Added the path=path keyword
;-

pro corr4flux, buf, x, y, data, nsig= nsig, reshape = reshape, $
			outer = outer, lemen=lemen, path=path
     
;; pfn = string('/3p/fuller/data/pc/'+buf+'.pc')
   pfn = concat_dir(path, buf+'.pc')			; 31-jan-94
   print, format = '(/,/," Reading file :", A, /)', pfn
   if (strlen(pfn) GT 19) then restore, file=pfn, /verbose $
   else return

;   a = fix(sqrt(float(N_ELEMENTS(data))));
   a = N_ELEMENTS(data);
   if N_ELEMENTS(reshape) then center_points, pcdata
   if N_ELEMENTS(outer) then outer_points, pcdata
;   pcdata = reform(float(pcdata.roi), a,a)
   time = pcdata.LIVE_TIME;
   pcdata = float(pcdata.roi);

   if N_ELEMENTS(lemen) then begin
	pcdata = fltarr(49);
	pcdata(0) = 633.0588;
	pcdata(1) = 907.9152;
	pcdata(2) = 1025.5130;
	pcdata(3) = 1022.5186;
	pcdata(4) = 1066.5046;
	pcdata(5) = 950.4421;
	pcdata(6) = 914.0438;
	pcdata(7) = 797.2570;
	pcdata(8) = 992.2537;
	pcdata(9) = 1046.2820;
	pcdata(10) = 1053.7740;
        pcdata(11) = 992.7240;
	pcdata(12) = 926.3719;
	pcdata(13) = 1060.5075;
	pcdata(14) = 882.8962;
	pcdata(15) = 969.8401;
	pcdata(16) = 1048.0777;
	pcdata(17) = 1036.4933;
        pcdata(18) = 1049.6776;
	pcdata(19) = 978.7311;
	pcdata(20) = 996.1217;
	pcdata(21) = 850.4812;
	pcdata(22) = 997.5146;
	pcdata(23) = 1086.9761;
	pcdata(24) = 973.0626;
	pcdata(25) = 1014.5046;
	pcdata(26) = 1042.6158;
	pcdata(27) = 981.4674;
	pcdata(28) = 650.2807;
	pcdata(29) = 831.6055;
	pcdata(30) = 1021.1853;
	pcdata(31) = 980.5006;
	pcdata(32) = 993.0566;
	pcdata(33) = 1003.9140;
	pcdata(34) = 845.3876;
	pcdata(35) = 398.0865;
	pcdata(36) = 639.8606;
	pcdata(37) = 781.6656;
	pcdata(38) = 836.5131;
	pcdata(39) = 793.6030;
	pcdata(40) = 769.6383;
	pcdata(41) = 702.4367;
	pcdata(42) = 218.8025;
	pcdata(43) = 376.0589;
	pcdata(44) = 536.8854;
	pcdata(45) = 657.6689;
	pcdata(46) = 620.3155;
	pcdata(47) = 528.4106;
	pcdata(48) = 410.1950;
    endif;

;
;  Get PC flux
;
   pcdata = pcdata/time;

;   XY = xy3(x(0:a-1,0), y(0,0:a-1))
   XY = fltarr(a, 2)
   XY(*,0) = x;
   XY(*,1) = y;

   sB = fltarr(6)	; Set up the coef array
   sB(0) = 1.0          ; Plateau Value
   sB(1) = 0.501	; Off-plateau slope
   sB(2) = 354.60	; Semi-major axis
   sB(3) = 304.86	; Semi-minor axis
;   sB(4) = 414.53       ; Centroid X
   sB(4) = 414.53+2.0   ; Centroid X -- corrected for Mons' report 6/20/93
   sB(5) = 423.35	; Centroid Y
;   sB(1) = 0.5		; Off-plateau slope
;   sB(2) = 355.0	; Semi-major axis
;   sB(3) = 310.0	; Semi-minor axis
;   sB(4) = 621.0       ; Centroid X
;   sB(5) = 415.0	; Centroid Y

   xx = x - sB(4);
   yy = y - sB(5);
   Re = sqrt((xx/sB(2))^2 + (yy/sB(3))^2)
   ie = where(Re LE 1.0, count)
   if count EQ 0 then avg = 1000. else avg = average(pcdata(ie))
   sB(0) = sB(0)*avg;
   sB(1) = sB(1)*avg;
   plateau2, XY, sB, sFit
;
;
   ii = fix(sqrt(N_ELEMENTS(sFit)));
   surface, reform(sFit-pcdata,ii,ii);
   if NOT(N_ELEMENTS(nsig)) then nsig = 2.0;
   err =  nsig*sqrt(float(pcdata))    ; 
   diff = abs( float(pcdata) - float(sFit))

   over = where(diff GE err);

   if (over(0) GE 0) then begin
	data(over) = data(over)*float(sFit(over))/float(pcdata(over));
	print, float(sFit(over))/float(pcdata(over));
   end;

   return
end;

;	(23-Mar-93)
;
;	Line Function -- To be called by CURVEFITA
;	Cone-Shaped Surface
;
PRO	LINE,R,A,F,PDER
;+
; NAME:
;	LINE
; PURPOSE:
;	Evaluate a function of a cone-shaped and optionally
;	return the value of its partial derivatives.
;	This function is used by CURVEFITA to fit this function and
;	a constant background to the data.
;
;	The routine CURVEFIT as supplied as an IDL user function has
;	been modified by Jim Lemen and subsequently by Rich Fuller
;       to compute a surface based on a singular dependent value, R,
;   	with an input of X and Y coordinates.
;
; CATEGORY:
;	E2 - SURFACE FITTING.
; CALLING SEQUENCE:
;	FUNCTA,X,A,F,PDER
; INPUTS:
;	A = PARAMETERS OF EQUATION DESCRIBED BELOW.
; OUTPUTS:
;	F = VALUE OF FUNCTION AT EACH R(I).
;
; OPTIONAL OUTPUT PARAMETERS:
;	PDER = (N_ELEMENTS(X),N_ELEMENTS(Y),4) ARRAY CONTAINING THE
;		PARTIAL DERIVATIVES.  P(I,J,K) = DERIVATIVE
;		AT ITH and JTH POINTS W/RESPECT TO KTH PARAMETER.
; COMMON BLOCKS:
;	NONE.
; SIDE EFFECTS:
;	NONE.
; RESTRICTIONS:
;	NONE.
; PROCEDURE:
;       F = A(2) + A(3)*R
;       R ===>  given (3/23/93, RAF)
; MODIFICATION HISTORY:
;       Written by R. Fuller, July 1992
;       Updated by R. Fuller, Jan 1993
;       Changed to input R directly by R Fuller, March 1993
;-
   R = R > (0.001);
   F = A(0) + A(1)*R

   PDER = FLTARR(N_ELEMENTS(R),2)

   PDER(0,0) = replicate(1.,n_elements(R))
   PDER(0,1) = R;

;   print, F, PDER

RETURN
END

;	(28-Mar-93)
;
;	Plateau Function -- To be called by CURVEFIT
;
PRO	PLATEAU2,XY,A,F,PDER
;+
; NAME:
;	PLATEAU
; PURPOSE:
;	This function is used by CURVEFIT to fit this function and
;	a constant background to the data.
;
; CATEGORY:
;	E2 - SURFACE FITTING.
; CALLING SEQUENCE:
;	PLATEAU,X,A,F,PDER
; INPUTS:
;	XY = 2-d array of values of independent variables = XY(N,2).
;	     The second argument is 0 for values of X and 1 for values of Y.
;	A = PARAMETERS OF EQUATION DESCRIBED BELOW.
; OUTPUTS:
;	F = VALUE OF FUNCTION AT EACH XY(I,K).
;
; OPTIONAL OUTPUT PARAMETERS:
;	PDER = (N_ELEMENTS(X),N_ELEMENTS(Y),6) ARRAY CONTAINING THE
;		PARTIAL DERIVATIVES.  P(I,J,K) = DERIVATIVE
;		AT ITH and JTH POINTS W/RESPECT TO KTH PARAMETER.
; COMMON BLOCKS:
;	NONE.
; SIDE EFFECTS:
;	NONE.
; RESTRICTIONS:
;	NONE.
; PROCEDURE:
;	X = XY(i,0) & Y = XY(i,1)
; MODIFICATION HISTORY:
;       Written by R. Fuller, August 1992
;       Modified by R. Fuller, March 1993 to take 2-d array of XY
;-

   X = XY(*,0) 				; Get X
   Y = XY(*,1)				; Get Y

   xx = x - A(4)
   yy = y - A(5)
   
   Re = sqrt((xx/A(2))^2 + (yy/A(3))^2)
   ie = where(Re LE 1.0, count)
   oe = where(Re GT 1.0)
   F = replicate(A(0), N_ELEMENTS(X));
   if (oe(0) GE 0) then F(oe) = A(0) + A(1)*(1.0-Re(oe))

   IF N_PARAMS(0) LE 3 THEN RETURN ;NEED PARTIAL?

   PDER = FLTARR(N_ELEMENTS(XY(*,0)), N_ELEMENTS(XY(*,1)),6)
   ptemp = fltarr(N_ELEMENTS(XY(*,0)), N_ELEMENTS(XY(*,1)));
   zz = replicate(0., N_ELEMENTS(XY(*,0)), N_ELEMENTS(XY(*,0)));

   PDER(0,0) = replicate(1., N_ELEMENTS(XY(*,0)), N_ELEMENTS(XY(*,0)));
   if count GT 0 then ptemp(ie) = zz(ie)
   if (oe(0) GE 0) then ptemp(oe) = -1.*Re(oe);
   PDER(0,1) = ptemp;
   if (oe(0) GE 0) then ptemp(oe) = A(1)*xx(oe)^2/A(2)^3/Re(oe);
   PDER(0,2) = ptemp;
   if (oe(0) GE 0) then ptemp(oe) = A(1)*yy(oe)^2/A(3)^3/Re(oe);
   PDER(0,3) = ptemp;
   if (oe(0) GE 0) then ptemp(oe) = A(1)*xx(oe)/A(2)^2/Re(oe)
   PDER(0,4) = ptemp;
   if (oe(0) GE 0) then ptemp(oe) = A(1)*yy(oe)/A(3)^2/Re(oe)
   PDER(0,5) = ptemp;

RETURN
END

pro plot_sxt_vig, high=high, low=low, ps=ps, path=path, noplot=noplot, qtest=qtest

ok = keyword_set(high) + keyword_set(low)
if ok ne 1 then begin
  message,'You must give the /high or the /low keyword',/cont
  return
endif

if n_elements(path) eq 0 then path = ''		; Default to current directory
set_plot,'X'					; For diagnostic plots

; -------------------------------
;  Check if any ws_* files exist
; -------------------------------
fff = findfile(concat_dir(path,'ws_*'),count=ncount)
if ncount eq 0 then begin
   print,'' & tbeep,5
   print,'** There are no wc_* files'
   print,'** Either you must specify the path= keyword,'
   print,'** or     you need to restore the data files from: ',findfile(concat_dir('$DIR_SXT_SENSITIVE','fuller_vignette.tar*'))
   print,'** Please refer to the documentation header'
   return
endif

; -------------------------------
;  Call the Fuller subroutines
; -------------------------------
if keyword_set(high) then begin
   eff_area_hi,ps=ps,path=path,noplot=noplot,qtest=qtest,fileout='cal_note37_hi.ps'
endif else begin
   eff_area_lo,ps=ps,/special,path=path,noplot=noplot,qtest=qtest,fileout='cal_note37_lo.ps'
endelse

end