pro eis_fix_pointing,file,hk_dir_head=hk_dir_head,fix_temp=fix_temp, $
       mhc_temp=mhc_temp,noaocs=noaocs
;+
;
;  Fixes fits files with erroneous pointins information....
;    ...assuming of course that one has access to the housekeeping files
;
;  $Id: eis_fix_pointing.pro 406 2010-07-04 11:33:51Z viggoh $
;
;-
  if n_params() ne 1 then begin 
    message,'eis_fix_pointing,file,hk_dir_head=hk_dir_head',/info
    return
  end
  aocs=0
  if n_elements(noaocs) eq 0 then aocs=1
  if n_elements(mhc_temp) eq 0 then mhc_temp=0
  if n_elements(fix_temp) eq 0 then fix_temp=0 else ccd_temp=1
;
  data=obj_new('eis_data',file)
  date_obs=data->getinfo('DATE_OBS') 
  date_end=data->getinfo('DATE_END')
  ccsds_last=data->getinfo('CCSDS_TI')
  nexp=data->getnexp()
; 
; possible bug here if rollover occurs during observation...
;  ...should be possible to catch on a day with a clear head....
;  ...which is not today.
;
  top_bits =ccsds_last 
  ti_1=*(data->getaux_data()).ti_1
  ti_2=*(data->getaux_data()).ti_2
  start_obs = string((top_bits/ulong(2)^28)*ulong(2)^28+ti_1[0]/16)
  end_obs = string((top_bits/ulong(2)^28)*ulong(2)^28+ti_2[n_elements(ti_1)-1]/16)
  t1ccsds=(top_bits/ulong(2)^28)*ulong(2)^28+ti_1/16.

  spcd=obj_new('eis_spcd',date_obs,date_end,start_obs,end_obs,dir_head=hk_dir_head)
  spcd->construct_pointing,aocs=aocs
  time_sc = ulong(reform(spcd->gettime_pt()))
  xcen_sc = spcd->getxcen()
  ycen_sc = spcd->getycen()
  sub=where(time_sc-shift(time_sc,-1) ne 0)
  if not spcd->getioerror() and n_elements(time_sc) ge 2 then begin
    time_sc=time_sc[sub]
    xcen_sc=xcen_sc[sub]
    ycen_sc=ycen_sc[sub]
    xcen = interpol(xcen_sc,time_sc,t1ccsds)
    ycen = interpol(ycen_sc,time_sc,t1ccsds)
  endif else begin
    xcen = fltarr(nexp)+xcen_sc[0]
    ycen = fltarr(nexp)+ycen_sc[0]
  endelse
  data->setxcen,xcen
  data->setycen,ycen
;
; subtract coarse mirror since eis_modfits uses the getxcen() function
;
;
  cal=obj_new('eis_cal',/quiet)
  cmirr_pos=data->getinfo('cmirr')
  fmirr=*(data->getaux_data()).fmirr
;
  crpix1 = 1
  crpix2 = 2
; pixel sizes
  cdelt1 = -1.*data->getdx_size() ; (arcsec pr pixel in solar-x)
  cdelt2 = 1.0 ; (arcsec pr pixel in solar-y)  
  crval1 = (spcd-> getcrval1()) + $
    (cal->getxmidmirrpos()-fmirr[0])*(cal->getfmirr_ss())*2.+ $
     0.5*data->getslit_width() + $
     (cmirr_pos-cal->getxmidcmirrpos())*(cal->getcmirr_ss()) + $
     (cal->getslit_pointing_offset(data->getslit_ind()))
  crval2 = (spcd-> getcrval2())-(cal->getymidslitpos())+(data->getys())[0]+ $
        (cal->getyoff_cmirr(data->date_obs2date(date_obs)))
;
  yw = data->getyw()
  nexp = data->getnexp()
  nexp_prp = data->getnexp_prp()
  nraster = 1
  if data->getinfo('fmir_ss') gt 0 then $
    nraster = nexp/nexp_prp ; nr of unique raster positions
  fovx = -1.*nraster*cdelt1
  fovy = yw*cdelt2
;
  xcen=data->getxcen()
;
  xpos=(spcd->getcrval1()) + $
    (cal->getxmidmirrpos()-fmirr[0])*(cal->getfmirr_ss())*2. + $
    (cmirr_pos-cal->getxmidcmirrpos())*(cal->getcmirr_ss()) + $
    (cal->getslit_pointing_offset(data->getslit_ind()))-fovx/2.
  ypos=(spcd->getcrval2()) + cdelt2*((yw+1)/2.-crpix2) - $
        (cal->getymidslitpos())+(data->getys())[0] + $
        (cal->getyoff_cmirr(data->date_obs2date(date_obs)))
  xcen0=xpos[0]
  ycen0=ypos[0]
  obj_destroy,cal
;
; UPDATE data header
;
  hdr = data->gethdr()
  main_hdr = hdr->getmain_hdr()
  bte_hdr = hdr-> getbte_hdr()
;
; main header values 
;
  fxaddpar,main_hdr,'XCEN',xcen0,'Center pointing (Solar X)'
  fxaddpar,bte_hdr,'XCEN',xcen0,'Center pointing (Solar X)'
  fxaddpar,main_hdr,'YCEN',ycen0,'Center pointing (Solar Y)'
  fxaddpar,bte_hdr,'YCEN',ycen0,'Center pointing (Solar Y)'  
  fxaddpar, main_hdr, 'CRPIX1', crpix1,'Reference pixel in solar-x'
  fxaddpar, main_hdr, 'CRPIX2', crpix2,'Reference pixel in solar-y'
  fxaddpar, main_hdr, 'CRVAL1', crval1[0],'Coordinate value of CRPIX1'  
  fxaddpar, main_hdr, 'CRVAL2', crval2[0],'Coordinate value of CRPIX2'  
  fxaddpar, main_hdr, 'CDELT1', cdelt1,'Pixel size in solar-x'
  fxaddpar, main_hdr, 'CDELT2', cdelt2,'Pixel size in solar-y' 
  fxaddpar, bte_hdr, 'CRPIX1', crpix1,'Reference pixel in solar-x'
  fxaddpar, bte_hdr, 'CRPIX2', crpix2,'Reference pixel in solar-y'
  fxaddpar, bte_hdr, 'CRVAL1', crval1[0],'Coordinate value of CRPIX1'  
  fxaddpar, bte_hdr, 'CRVAL2', crval2[0],'Coordinate value of CRPIX2'  
  fxaddpar, bte_hdr, 'CDELT1', cdelt1,'Pixel size in solar-x'
  fxaddpar, bte_hdr, 'CDELT2', cdelt2,'Pixel size in solar-y'
;
; max and min values of pointing has changed.
; TDMIN and TDMAX must be set again
; First find which column the pointings are stored in.
;
  naxis1=long(fxpar(bte_hdr,'NAXIS1'))
  for i=1,50 do begin
    nr=strtrim(string(i),2)
    param=string(data->getinfo('TTYPE'+nr))
    case param of 
    'XCEN_TI1': begin
       xcen_nr=nr
       form=strtrim(fxpar(bte_hdr,'TFORM'+nr),2)
       if strmid(form,strlen(form)-1) eq 'D' then begin
         param1='TFORM'+nr
         fxaddpar,bte_hdr,param1,'1E','Real*4 (floating point)'
         naxis1=naxis1-4l
       endif
               end
    'YCEN_TI1': begin
       ycen_nr=nr
       form=strtrim(fxpar(bte_hdr,'TFORM'+nr),2)
       if strmid(form,strlen(form)-1) eq 'D' then begin
         param1='TFORM'+nr
         fxaddpar,bte_hdr,param1,'1E','Real*4 (floating point)'
         naxis1=naxis1-4l
       endif
               end
    else:
    endcase
  endfor
;
  param1 = 'TDMIN'+xcen_nr
  param2 = 'TDMAX'+xcen_nr
  fxaddpar, bte_hdr, param1, min(xcen)
  fxaddpar, bte_hdr, param2, max(xcen)
  param1 = 'TDMIN'+ycen_nr
  param2 = 'TDMAX'+ycen_nr
  fxaddpar, bte_hdr, param1, min(ycen)
  fxaddpar, bte_hdr, param2, max(ycen)
;
  if fix_temp then begin
    spcd->construct_ccdtemp
    time_sc = double(reform(spcd->gettime_sts1()))
    ccda_temp = double(spcd->getccda_temp())
    ccdb_temp = double(spcd->getccdb_temp())
    if not spcd->getioerror() and n_elements(time_sc) ge 2 then begin
      ccda_temp = float(interpol(ccda_temp,time_sc,double(t1ccsds)))
      ccdb_temp = float(interpol(ccdb_temp,time_sc,double(t1ccsds)))
      ccda_temp = eis_convert_ccda_temperature(ccda_temp)
      ccdb_temp = eis_convert_ccdb_temperature(ccdb_temp)
    endif else begin
      if not spcd->getioerror() then begin
        ccda_temp = eis_convert_ccda_temperature(ccda_temp[0])+fltarr(nexp)
        ccdb_temp = eis_convert_ccda_temperature(ccdb_temp[0])+fltarr(nexp)
      endif else begin
        ccda_temp = fltarr(nexp)
        ccdb_temp = fltarr(nexp)
      endelse
    endelse
    data->setccda_temp,ccda_temp
    data->setccdb_temp,ccdb_temp
    for i=1,50 do begin
      nr=strtrim(string(i),2)
      param=string(data->getinfo('TTYPE'+nr))
      case param of 
      'CCDAT_TI': begin
         ccdat_nr=nr
         form=strtrim(fxpar(bte_hdr,'TFORM'+nr),2)
         if strmid(form,strlen(form)-1) eq 'D' then begin
           param1='TFORM'+nr
           fxaddpar,bte_hdr,param1,'1E','Real*4 (floating point)'
           naxis1=naxis1-4l
         endif
                 end
      'CCDBT_TI': begin
         ccdbt_nr=nr
         form=strtrim(fxpar(bte_hdr,'TFORM'+nr),2)
         if strmid(form,strlen(form)-1) eq 'D' then begin
           param1='TFORM'+nr
           fxaddpar,bte_hdr,param1,'1E','Real*4 (floating point)'
           naxis1=naxis1-4l
         endif
                 end
      else:
      endcase
    endfor
    param1 = 'TDMIN'+ccdat_nr
    param2 = 'TDMAX'+ccdat_nr
    fxaddpar, bte_hdr, param1, min(ccda_temp) 
    fxaddpar, bte_hdr, param2, max(ccda_temp)
    param1 = 'TDMIN'+ccdbt_nr
    param2 = 'TDMAX'+ccdbt_nr
    fxaddpar, bte_hdr, param1, min(ccdb_temp)
    fxaddpar, bte_hdr, param2, max(ccdb_temp)
;
    if mhc_temp then begin
      spcd->construct_mhctemp
      s=spcd->getmhc_temp()
      t15=double(eis_convert_mhc_temperature(s.mhc_hz_t15,s.mhc_ref_ther_0,s.mhc_cal_ther_0))
      ti5=smooth(t15,11<(size(t15))[1],/edge)
      t10=double(eis_convert_mhc_temperature(s.mhc_hz_t10,s.mhc_ref_ther_0,s.mhc_cal_ther_0))
      ti0=smooth(t10,11<(size(t10))[1],/edge)
      if not spcd->getioerror() and n_elements(s.time_sc) ge 2 then begin
        t10 = float(interpol(double(t10),double(s.time_sc),double(t1ccsds)))
        t15 = float(interpol(double(t15),double(s.time_sc),double(t1ccsds)))
      endif else begin
        if not spcd->getioerror() then begin
          t10 = float(t15[0])+fltarr(nexp)
          t15 = float(t15[0])+fltarr(nexp)
        endif else begin
          t10 = fltarr(nexp)
          t15 = fltarr(nexp)
        endelse
      endelse
      data->setmhc_hz_t10,t10
      data->setmhc_hz_t15,t15
;
      mhc_hz_t10nr=-1
      for i=1,50 do begin
        nr=strtrim(string(i),2)
        param=string(data->getinfo('TTYPE'+nr))
        case param of 
        'MHC_HZ_T10': begin
           mhc_hz_t10nr=nr
                   end
        'MHC_HZ_T15': begin
           mhc_hz_t15nr=nr
                   end
        else:
        endcase
      endfor
;
;  ADD mhc_hz_t10 and mhc_hz_15 column if not already there
;
      if mhc_hz_t10nr eq -1 then begin       
        fxbaddcol, col, bte_hdr, float(t10[0]), 'MHC_HZ_T10', $
                  'MHC HZ TEMP 10 (MHC_T10(TI) )', $
                  tunit = 'C',$
                  tdmin=min(t10), $
                  tdmax=max(t10)
        mhc_hz_nr=ccdbt_nr+1
        naxis1=naxis1+4l      
        fxbaddcol, col, bte_hdr, float(t15[0]), 'MHC_HZ_T15', $
                  'MHC HZ TEMP 15 (MHC_T15(TI) )', $
                  tunit = 'C',$
                  tdmin=min(t15), $
                  tdmax=max(t15)
        mhc_hz_t10nr=ccdbt_nr+1
        mhc_hz_t15nr=ccdbt_nr+2
        naxis1=naxis1+4l
      endif else begin
        param1 = 'TDMIN'+mhc_hz_t10nr
        param2 = 'TDMAX'+mhc_hz_t10nr
        fxaddpar, bte_hdr, param1, min(t10) 
        fxaddpar, bte_hdr, param2, max(t10)
        param1 = 'TDMIN'+mhc_hz_t15nr
        param2 = 'TDMAX'+mhc_hz_t15nr
        fxaddpar, bte_hdr, param1, min(t15) 
        fxaddpar, bte_hdr, param2, max(t15)
      endelse
    endif
  endif
;  
  fxaddpar,main_hdr,'NAXIS1',naxis1,'Number of bytes per row'
  fxaddpar,bte_hdr,'NAXIS1',naxis1,'Number of bytes per row'
  hdr-> setmain_hdr, main_hdr
  hdr-> setbte_hdr, bte_hdr
;
; SAVE data to fits file
;
  newfile=data->getfilename()
;
; take '.gz' off name if there is one
;
  newfile=strsplit(newfile,'\.gz',/extract,/regex)
  data->save,file=newfile
;
return
end