Pro sedioc

;NAME:sedioc.pro
;
;PURPOSE:This program is to find the best CSMM (Scan Mirror) position for 
;        SED mode and find the tilt between the array and the slit.
;
;CALLING SEQUENCE:.r sedioc.pro
;                 sedioc
; Will call tilt.pro
;
;INPUTS:name of file with a list of the fits filenames. The first line of this
; file should be the number of names.
;  The fits files should be in the SED-IOC format: 
;              background, stim 12 on/off pairs, stim 
;  Each fits file is a telescope position. Each extension is a stim, background
; or new CSMM(scan mirror) positon
; 
;OPTIONAL KEYWORD:none
;
;OUTPUTS: seddetails - a file with the max col for each CSMM pos, the max 
; csmm pos, the max telescope position, etc.
; CSMMfit - details about fitting the curve
;
;SIDE EFFECTS:will change the color scheme to Rainbow+white
;
;RESTRICTIONS:
;
;EXAMPLES:
;
;MODIFICATION HISTORY: 1/8/03 added headers to the sedIOC files so the 
;plots can also have the scan mirror position in counts
;
;------------------------------------
Print, ''
PRINT, 'Remember to change from counts to arcsec in real data!'
Print, ''
units:
sm_units='counts'
Print, 'Should the scan mirror units be in counts or arcsecs? c/a'
read, sm_units
if (sm_units ne 'c' and sm_units ne 'a') then begin
  print, 'Type c or a to choose counts or arcsecs for the scan mirror units'
  goto, units
endif
if (sm_units eq 'c') then sm_units='counts'
if (sm_units eq 'a') then sm_units='arcsec'

;set the color right
loadct,39

;Each file should be a new telescope position, each fits exention (pic) 
;should be a new Scan Mirror (CSMM) position
filename=''
PRINT, 'First line of the file should be the number of files in the list.'
READ, filename, PROMPT='Enter the name of the file with the list of names:'
pmax=1  ;number of files in filename
tempname=''
Close, 5
OPENR, 5, filename 
READF,5,pmax
maxtele=fltarr(pmax,2)
name=STRARR(pmax)
for p=0, pmax-1 do begin
  READF,5, tempname
  name(p)=tempname
endfor 
Close, 5



;File that will hold all the detailed info
CLOSE, 10
OPENW, 10, 'seddetails'
apercol=fltarr(pmax,12)   ;sum over just that aperture for each tel. pos
again=''
smpos=fltarr(pmax,12)
smMax=fltarr(pmax)
RAS=fltarr(pmax)
DECS=fltarr(pmax)

;File with fitting info
CLOSE, 20
OPENW,20, 'CSMMfit'

;file with tilt info
CLOSE,15
Openw,15, 'Tilt_details'
PRINTF,15, 'Tilt_details'


For p=0,pmax-1  do begin ;each telescope position(fits file)

;make sure the number of cycles is correct
;ncycles=FLOAT(hvalue(name(p),'NCYCLES',EXTEN=0))
;if (ncycles ne 12.0) then begin
;PRINT, 'Number of cycles not valid'
;STOP
;ENDIF

;Read from the header what the telescope position is.
;The RA's will be in array RAS and the Declination will be in DECS
;RAS(p)=FLOAT(hvalue(name(p),'CRVAL1',Exten=0))
;DECS(p)=FLOAT(hvalue(name(p),'CRVAL2',Exten=0))
;both are in degrees


PRINT,''
PRINT, 'On file: ',name(p)
cube=rdfitscube(name(p))


PRINTF, 10, ' '
PRINTF, 10, 'For file:  ', name(p)
PRINTF, 10, ' '
PRINTF, 10, '     DCE#,  Highest value,   col#,    sum over aperture'
PRINTF, 10, '               DN/sec                      DN/sec      '


;Subtract off the dark sky frames 
;First 2 extensions(planes) of cube are a background and a stim
;each pic is a CSMM position
pic=fltarr(32,32,12)
for i=0,11 do begin
   pic(*,*,i)= cube(*,*,(2*i+2))-cube(*,*,(2*i+3))
endfor

;to get the aperture
;get the median of all 12 pics
sumpic=fltarr(32,32)

for i=0,31 do begin
 for j=0,31 do sumpic(i,j)=MEDIAN(pic(I,j,*),/EVEN)
endfor

sumpiccol=TOTAL(sumpic, 2)
Plot, sumpiccol , TITLE='Median of all Pictures',$
     ytit='Median of Fluxes in DN/sec', $
     xtit='Column', background=255, color=0,xticks=8, $
      xticklen=1, xgridstyle=1,xrange=[0,32],xstyle=1

;ask for input
ap:
aper1=1
aper2=1
READ, aper1, PROMPT='At what column should the aperture start?'
READ, aper2, PROMPT='At what column should the aperture stop?'
if (aper1 lt 0 or aper2 gt 23) then begin
print, "The range of the aperture should be between 0 and 23."
goto, ap
endif

;for the 12 on/off pairs
cols=fltarr(12,32) ; the sum of the col,for each pic
maxcol=fltarr(12)
bkgrd=fltarr(12,32)  ; background for that col
bxs=fltarr(23-aper2+aper1)
bys=fltarr(23-aper2+aper1)
bxys=fltarr(23-aper2+aper1)
bxsqr=fltarr(23-aper2+aper1)

;1Set up the plot
plot,[0], /nodata,xrange=[0,33], yrange=[200,150000], $ 
   tit='Sums of Columns for ' +name(p), xtit='column #',ytit='DN/sec', $
   background=255, color=0,xticklen=1, xgridstyle=1,xticks=8
pcolor=[20,40,70,100,150,190,20,40,70,100,150,190 ]

;for each CSMM position for telescope position p
for picnum=0,11 do begin
  cols(picnum,*)=TOTAL(pic(*,*,picnum),2)


;getting the background for each col:
jj=0
  for j=0,(aper1-1) do begin
     bkgrd(picnum,jj)=cols(picnum,j)
     bxs(jj)=j
     jj=jj+1
  endfor

  for j=(aper2+1),23 do begin
      bkgrd(picnum,jj)=cols(picnum,j)
      bxs(jj)=j
      jj=jj+1
  endfor

  for j=0,jj-1 do begin
  bys(j)=bkgrd(picnum, j)
  bxys(j)=bxs(j)*bys(j)
  bxsqr(j)=bxs(j)*bxs(j)
  endfor

delta=(jj*TOTAL(bxsqr)-(TOTAL(bxs))^2)
;slope
bm=(jj*TOTAL(bxys)-TOTAL(bxs)*TOTAL(bys))/delta
;y-int
bb=(TOTAL(bxsqr)*TOTAL(bys)-TOTAL(bxs)*TOTAL(bxys))/delta

for j=0,31  do bkgrd(picnum,j)=bm*j+bb
;subtracting off background
;If (p eq 1) then begin
;PRINT, 'checking background'
;STOP
  cols(picnum,*)=cols(picnum,*)-bkgrd(picnum,*)
;endif else begin
;  cols(picnum,*)=cols(picnum,*)-bkgrd(picnum,*)
;endelse

;get the scan mirror position in counts (?)
switch_units=''
count_start:
if (sm_units eq 'counts') then begin
 test_SM=FLOAT(hvalue(name(p),'SM_POS ',EXTEN=1))
 if (test_SM eq -999) then begin 
 PRINT, ' '
 PRINT, 'The postion in counts cannot be found in the header.'
 READ,switch_units, PROMPT='Would you like to switch to arcsec? y/n: '
   if (switch_units eq 'y') then begin
       sm_units='arcsec'
       GOTO, arcsec_start
   endif
 endif

  ii=0
  FOR i=3,25,2 do begin ;only the on positions
   smpos(p,ii)=FLOAT(hvalue(name(p),'SM_POS ',EXTEN=i))
   ii=ii+1
  endfor
endif else begin
arcsec_start: 
;get scan mirror postion in arcsec
test_SM=FLOAT(hvalue(name(p),'CSM_POS ',EXTEN=1))
 if (test_SM eq -999) then begin
 PRINT, ' '
 PRINT, 'The postion in arcsecs cannot be found in the header.' 
 READ,switch_units, PROMPT='Would you like to switch to counts? y/n: ' 
   if (switch_units eq 'y') then begin
        sm_units='counts'
        GOTO, count_start
   endif
 endif

 ii=0
  FOR i=3,25,2 do begin ;only the on positions
   ;Real name, will be in arcsec
   smpos(p,ii)=FLOAT(hvalue(name(p),'CSM_POS ',EXTEN=i))
   ii=ii+1
  endfor
endelse


; Ploting
if picnum lt 6 THEN ls=2 else ls=3
  OPLOT, cols(picnum,*), color=pcolor[picnum], linestyle=ls

;Labels
;yleg=10^(0.98^((picnum+1)/1.1)*!y.crange[1])
yleg=(0.98^((picnum+1)/1.1)*!y.crange[1])
xyouts, 20,yleg, 'CSMM Position'+STRING(picnum+1),  color=0

oplot, [30,35], [yleg,yleg], color=pcolor[picnum],linestyle=ls


;adding just the aperture
for i=aper1,aper2 do apercol(p,picnum)=apercol(p,picnum)+cols(picnum,i)

  maxcol(picnum)=MAX(cols(picnum,*), m1)
  PRINTF, 10,  picnum+1, maxcol(picnum), m1,apercol(p,picnum)
endfor  ;for picnum

PRINT, 'Sanity Check!'
PRINT, 'Type .c to continue'
STOP     

;To find the best CSMM Position

PRINT,''
PRINT,'Finding best CSMM position'
PRINT,''

;start with fitting all positions
fit1=0
fit2=11
PosLoop:

;Fit gaussian to apercol(p,*)

x=DINDGEN(fit2-fit1+1)
for j=0,(fit2-fit1) do x(j)=fit1+j
smpos2=fltarr(fit2-fit1+1)
  for j=0,fit2-fit1 do smpos2(j)=smpos(p,x(j))
A=fltarr(4)
aperfit=MPFITPEAK(smpos2,apercol(p,fit1:fit2),A)
x2=fltarr(n_elements(smpos2)*4)
x2(0)=smpos2(0)
xrange=smpos2(n_elements(smpos2)-1)-smpos2(0)
for j=1,n_elements(x2)-2 do x2(j)=x2(0)+j*(xrange/n_elements(x2))
x2(n_elements(x2)-1)=smpos2(n_elements(smpos2)-1)
;z=(smpos2-A(1))/A(2)
z=(x2-A(1))/A(2)


PRINT,''
PRINT, 'Coeffs for original: ',A
PRINT,''

smMax(p)=A(1)
apercheck=A(0)*EXP(-(z^2)/2)+A(3)

Plot, smpos2, apercol(p,fit1:fit2),xtit='CSMM Position in '+sm_units,$
    ytit='Flux thru Aperture in DN/sec',background=255,color=0,$
    tit='Fit and Original for CSMM Position'
OPLOT, x2,apercheck, linestyle=2,color=70
;LABELS
yleg=(0.98^(1/1.1)*!y.crange[1])
xleg=!x.crange[0]+(!x.crange[1]-!x.crange[0])*.75
xyouts, xleg,yleg,'Original Data',color=0
yleg=(0.98^(3/1.1)*!y.crange[1])
xyouts,xleg,yleg,'Fit of Original',color=70


maxcsmmpos=MAX(apercol(p,*), csmm)
PRINT, 'The best (theoretical) CSMM pos is; ',smMax(p)
Print, 'The CSMM position with the highest flux is: '$
   +'# '+strtrim(String(csmm+1),2)+', '+strtrim(string(smpos(p,csmm)),2)+$
   '  counts'
Print, 'The highest flux is:', maxcsmmpos,' DN/sec'

READ, again, PROMPT='Is this fit ok? y/n:'
If (again eq 'n') Then begin
 oops:
   Plot, apercol(p,*), xtit='CSMM Position #', $
     ytit='Sum over aperture in DN/sec',$
     tit='Flux through aperture vs. CSMM Position #',background=255, color=0 $
     ,xticklen=1, xgridstyle=1,xticks=8

   PRINT, 'Choose which columns to be fit (>3) ,'
   READ, fit1,PROMPT='Start fit at what position?'
   READ, fit2, PROMPT='End fit at what position?'
     if (fit2-fit1 le 3) then begin
     print, 'There must be more than 3 columns to fit.'
     goto, oops
     endif

   Goto, PosLoop
ENDIF


;to find column with most flux
;theoretical:

PRINT,''
PRINT, 'Finding the tilt of the array '
PRINT, 'and where the object is along the slit.' 
PRINT,''
;READ, goodCSMM, PROMPT='Pick a CSMM position with good flux:  '
;picks col with highest flux
goodCSMM=csmm

ColLoop:
CLOSE,15
OPENW,15, 'Tilt_details',/APPEND
PRINTf,15,' '
PRINTF,15, 'For file: ',name(p)
Close,15
;calling tilt.pro to get the angle of the tilt
slitpos=tilt(pic(*,*,csmm))
goodmax=slitpos(1)

;absolute:
maxspc=MAX(sumpiccol, mspc)
PRINT, 'The highest flux is on column:', mspc

READ, again, PROMPT='Would you like to fit the tilt again? y/n:'
;Plot, apercol(p,*), xtit='CSMM Position #', $
;   ytit='Sum over aperture in DN/sec',$
;   tit='Flux through aperture vs. CSMM Position #',background=255, color=0 $
;   ,xticklen=1, xgridstyle=1,xticks=8
IF (again eq 'y') THEN GOTO, ColLoop

PRINTF, 10, '      Max Col,   Observed Max'
PRINTF, 10, goodMax, mspc


PRINTF, 20, 'For file', name(p)
PRINTF, 20, 'Fit for CSMM position: '
PRINTF, 20, 'Coefficients:    ', A
PRINTF, 20, 'CSMM Max:     ',smMax(p), csmm,smpos(p,csmm)
PRINTF, 20, ''
PRINTF, 20, 'Pos of Slit: ', goodMax, mspc
PRINTF, 20, ''
PRINTF, 20, ''

PRINTF,10, ''
PRINTF,10, 'The CSMM position with the highest flux is:',$
           csmm,smpos(p,csmm),'counts'
PRINTF,10, 'The highest flux (DN/sec) is:', maxcsmmpos
PRINTF,10, 'The calculated best CSMM position is: ',smMax(p),' counts'

PRINTF,10 , ' '
PRINTF,10, 'Tilt info:'
PRINTF,10, 'Slope: ',slitpos(0),'  +/-',slitpos(3)
PRINTF,10, 'Y-Int= ',slitpos(1),'  +/-',slitpos(4)
PRINTF,10, 'Angle= ',slitpos(2)
PRINTF,10 , ' '
;Compare the movement of the telescope to the movement of the scan mirror
;For i=0,pmax-2
;deltasm=smMax(i)-smMax(i+1)
;GCIRC,1,RAS(i)/15.0,DECS(i),RAS(i+1)/15.0,DECS(i+1),deltatel
;;deltatel is in arcsec
;;RA's should be in hours, so they are divided by 15
;PRINT, 'Between positon',i,' and ',i+1
;Print,'Change in scan mirror is', deltasm
;Print,'Change in telescope position is', deltatel
;
;PRINTF,10, 'Between positon',i,' and ',i+1
;Printf,10,'Change in scan mirror is', deltasm
;Printf,10,'Change in telescope position is', deltatel
;Endfor
Print, 'This will tell the angle that the scan mirror moved and'
PRINT , 'the angle that the telescope moved between each positon.'

maxtele(p)=smMax(p)

ENDFOR   ; p - every fits file

mtele=MAX(maxtele, mt)

;PRINTF,10, ' '
;PRINTF,10, 'Max telescope position is ', mt+1
;PRINT, 'Max telescope position:',mt+1
 CLOSE, 10
 CLOSE,15
 CLOSE, 20
 PRINT, 'Details have been printed to file seddetails.'



    s=Max(smpos,pos2,MIN=pos1)
     xleg=pos1+(s-pos1)*.7
     xleg1=pos1+(s-pos1)*.9
     xleg2=pos1+(s-pos1)*.99
PLOT, [0], /nodata, xrange=[pos1,s],yrange=[2000,maxcsmmpos], $
  xtit='CSMM Positions in Counts',ytit='Sum Over Aperture in DN/sec',$
  tit='Best CSMM Postions',background=255,color=0

 for p=0, pmax-1 do begin
    OPLOT, smpos(p,*),apercol(p,*),color=pcolor(p), linestyle=p
   ; yleg=(0.98^((p+1)/1.1)*!y.crange[1])
   ; yleg=(0.98^((p+1)/1.1)*maxcsmmpos)
     yleg=maxcsmmpos/1.5 - p*.05*maxcsmmpos
    xyouts, xleg,yleg, 'Telescope Position '+strtrim(STRING(p+1),2),$
     color=pcolor(p)
    oplot, [xleg1,xleg2], [yleg,yleg],color=pcolor(p),linestyle=p
  endfor

    



end