;+
; NAME:
;     OVSA2MAP
; PURPOSE:
;     Routine to read an OVSA FITS image file and return a
;     MAP object with correct rotation for solar P angle.
; CATEGORY:
;     OVRO MAPPING
; CALLING SEQUENCE:
;     map = ovsa2map(filename[,header])
; INPUTS:
;     filename  A string containing path and name of .FTS
;                 file containing the OVSA image.
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
; ROUTINES CALLED:
;     sxpar, fits2map, rot_map
; OUTPUTS:
;     map       The map object containing the image and
;                 associated information about it.
;     header    Optional output containing a verbatim copy
;                 of the FITS file header
; COMMENTS:
; SIDE EFFECTS:
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 28-Feb-2003 by Dale E. Gary
;     28-Apr-2009  DG
;       Added scaling by TBMAX
;     27-Jan-2014 GN
;       Added code to return correct integration time as duration
;     05-Feb-2014 GN
;       Added code to return standard fits2map output if not an OVSA fits file
;     07-May-2014 GN
;       Added offset tag to ovsa map structure, 
;       which is intended to indicate that a global shift of the map center may be needed
;     09-May-2014 GN
;       Added beam prmeter tags to the ovsa map structure
;     11-Oct-2014  DG
;        Added a 1/2-pixel offset to the southwest because the definition of
;        the map center should be the center of a pixel.  The correct offset
;        was determined by "imaging" the beam and observing the offset as a
;        function of frequency.
;        ;-

function ovsa2map,file,header,freq_array=freq_array,offset=offset, no_offset=no_offset
   fits2map,file,ovsa,head=header
   if strcompress(fxpar(header,'INSTRUME'),/rem) eq 'RHESSI' then hsi_fits2map,file,ovsa
   instr=fxpar(header,'INSTRUME')
   if size(instr,/tname) eq 'STRING' then begin
   if instr eq 'OVSA' then begin
      ovsa.dur=anytim(sxpar(header,'TEND'))-anytim(sxpar(header,'TSTART'))
      if n_elements(freq_array) eq 0 then freq_array=[1.20000, 1.40000, 1.60000, 1.80000, 2.00000, 2.40000, 2.60000, 2.80000, $
      3.20000, 3.40000, 3.60000, 3.80000, 4.20000, 4.40000, 4.60000, 4.80000, 5.00000, $
      5.20000, 5.40000, 5.60000, 6.20000, 6.60000, 7.00000, 7.40000, 7.80000, 8.20000, $
      8.60000, 9.00000, 9.40000, 10.0000, 10.6000, 11.2000, 11.8000, 12.4000, 13.2000, $
      14.0000, 14.8000, 15.6000, 16.4000, 18.0000]
       bif = sxpar(header,'BIF')
       eif = sxpar(header,'EIF')
       poln = sxpar(header,'POLN')
       freq=minmax(freq_array[bif:eif])
       if freq[0] ne freq[1] then $
       freq=strcompress(string(freq,format="(f4.1,'-',f4.1,'GHz')"),/rem) else $
       freq=strcompress(string(freq[0],format="(f4.1,'GHz')"),/rem)
       if ovsa.dur ne 0 then dur=strcompress(string(ovsa.dur,format="(i10,'s')"),/rem) else dur=''
       ovsa.id=ovsa.id+' '+poln+' '+freq+' '+dur
       tbmax = fxpar(header,'TBMAX')
  	   if (tbmax ne 0) then ovsa.data = ovsa.data*tbmax/max(ovsa.data)
  	   solp = sxpar(header,'SOLP')
  	   ovsa.roll_angle = -solp
  	   ovsa.roll_center = [0,0]
  	   ovsa = rot_map(ovsa,solp,center=[0,0])
  	  
  	   ; Read info from the header to retrieve the 1/2 power beamsize
  	   a = sxpar(header,'BMAJ')*sqrt(alog(2.))  ; Semi-major half-power axis
  	   b = sxpar(header,'BMIN')*sqrt(alog(2.))  ; Semi-minor half-power axis
  	   phi = sxpar(header,'BPOS')               ; Position angle of major axis
  	   solp = sxpar(header,'SOLP')              ; Solar P angle
  	   ; Correct angle for solar P angle (since image is rotated by P angle)
  	   phi = (phi - solp)*!dtor
  	   ovsa=create_struct(ovsa,'a',a,'b',b,'phi',phi)
  
       ; Shift map center to SW by 1/2 pixel in x and y
       ; to account for definition of center of map as center
       ; of the central pixel.
       ovsa.xc -= ovsa.dx/2.
       ovsa.yc -= ovsa.dy/2.
       
  	   test_offset = sxpar(header,'offset')
  	   if n_elements(test_offset) lt 2 then begin
    	   if n_elements(offset) ne 2 then offset=[0.0,0.0]
    	     offset=float(offset)
      	   if offset[0] eq 0 and offset[1] eq 0 then off=0b else off=1b
      	   ovsa=create_struct(ovsa,'offset',offset,'off',off)
      	   ovsa.xc+=ovsa.offset[0]
      	   ovsa.yc+=ovsa.offset[1] 
      	  if keyword_set(no_offset) then begin
      	    ovsa.xc=ovsa.xc-ovsa.offset[0]
            ovsa.yc=ovsa.yc-ovsa.offset[1]
            ovsa.off=0b
      	  endif
    	 end  
    	 
	   end
	 end
   return,ovsa
end