; file: makedopsumcubes.pro for 2007-04-12/QuietSun/Cubes
; last: Oct 15 2010 
; note: assocs into 2 line-wing cubes to produce Dopplermap and sum intensity

; input file parameters 
inpath='/home/rutten/data/DOT/2007-04-12/QuietSun/Cubes/'
inheader='QuietSun-'
intail='-20070412-cube-1034x776x180-dt30'

; output file parameters 
outpath='/home/rutten/data/DOT/2007-04-12/QuietSun/Cubes/'
outpath='./'
outheader=inheader
outtail=intail

; wavelength selection 
inwavident=['ha-0500','ha+0500']
outwavident=['ha-dopp','ha-wing']

; input cube dimensions
nx=1034
ny=776
nt=180
dt=30     ; cadence in seconds

; output cube dimensions = the same
nx2=nx
ny2=ny
nt2=nt
dt2=dt


; ================ this program should be parameter-free from here

; open and assoc cube1
get_lun, unit1
openr,unit1,inpath+inheader+inwavident[0]+intail
p1 = assoc(unit1, intarr(nx,ny))
  
; open abd assoc cube2
get_lun, unit2
openr,unit2,inpath+inheader+inwavident[1]+intail
p2 = assoc(unit2, intarr(nx,ny))

; goto,SKIPDOP   ; insert to skip production of the Dopplergram  cube

; make Doppler cube
cube=fltarr(nx,ny,nt)
for it=0,nt-1 do begin
  a=float(p1[it])
  b=float(p2[it])
  cube[*,*,it]=(a-b)/(a+b)
endfor
cube=histo_opt(cube,0.001,0.001)    ; Suetterlin pro to take out extremes 

; scale cube to within integer wordlength while maintaining zero
maxdop=max(cube)
mindop=min(cube)
print,' mindop, maxdop =',mindop,maxdop
if (-mindop gt maxdop) then maxdop=-mindop
cube=fix(cube*30000./maxdop)

; write Doppler cube
get_lun, unit
openw,unit,outpath+outheader+outwavident[0]+outtail
writeu,unit,cube
free_lun,unit
print,'  wrote ',outheader+outwavident[0]+outtail

SKIPDOP:

; make wing sum cube
cube=fltarr(nx,ny,nt)
for it=0,nt-1 do cube[*,*,it]=float(p1[it])+float(p2[it])
cube=histo_opt(cube,0.001,0.001)  ; remove extremes
maxcube=max(cube)

; scale cube to within integer wordlength
cube=fix(cube*30000./maxcube)
  
; write wing sum cube
get_lun, unit
openw,unit,outpath+outheader+outwavident[1]+outtail
writeu,unit,cube
free_lun,unit
print,'  wrote ',outheader+outwavident[1]+outtail

; close the input files
free_lun,unit1
free_lun,unit2

end