;+
; NAME:
;    wimagr_fits2map
; PURPOSE:
;    OVSA INDEX to MAP conversion routine
; CATEGORY:
;     OVSA
; CALLING SEQUENCE:
;     wimagr_index2map,data,index,info
; INPUTS:
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
;
; ROUTINES CALLED:
;
; OUTPUTS:
; COMMENTS:
; SIDE EFFECTS:
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 15-May-2014 Gelu M. Nita
;     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 wimagr_index2map,data,index,info
  nmaps=n_elements(index)
  nstrings=n_elements(index[0].header)
  headers=reform(index.header,nstrings,info.npol,info.nfreq,info.ntim)
  for i=0,info.npol-1 do  begin
    for j=0,info.nfreq-1 do  begin
      for k=0,info.ntim-1 do  begin
        header=headers[*,i,j,k]
        index2map,fitshead2struct(header),data[*,*,i,j,k],ovsa
        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 display 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
        offset=float(info.offset)
        ovsa.xc+=offset[0]
        ovsa.yc+=offset[1]

        ; 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.

        if offset[0] eq 0 and offset[1] eq 0 then off=0b else off=1b
        ovsa=create_struct(ovsa,'a',a,'b',b,'phi',phi,'offset',offset,'off',off)
        ovsa.dur=anytim(sxpar(header,'TEND'))-anytim(sxpar(header,'TSTART'))
        if tag_exist(info,'freq') then freq_array=info.freq $
        else 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)
        if i+j+k eq 0 then map_arr=ovsa else map_arr=[map_arr,ovsa]
      endfor
    endfor  
   endfor
    return,map_arr
  end