;
;
;
;
;
;
;
;

pro dem, files=files

IF NOT keyword_set(files) THEN BEGIN
    files=''
    read,prompt='Enter AIA image identifier (i.e. "master"):  ',files
ENDIF

;get AIA response functions
aia_resp=aia_get_response(/dn)
filters=['A131','A171','A193','A211','A335','A94']

;set exposure times, identically 1s based on pre-processing of de-rotation
exp_times=fltarr(6)+1

;set directory pathing information and find master de-rotated files
path='/home/aldebaran/schonfsj/research/f10.7/'
master_filenames=strcompress(path+'*'+string(files)+'.fits',/remove_all)
master_files=file_search(master_filenames)
;remove 304 A band for the moment
master_files=master_files(where(strpos(master_files,'304') eq -1))
print,master_files

;define data array for the de-rotated images
master=fltarr(4096,4096,n_elements(master_files))
;master=fltarr(512,512,n_elements(master_files))

;populate data array
FOR i=0,n_elements(master_files)-1 DO master[*,*,i]=mrdfits(master_files[i])
;FOR i=0,n_elements(master_files)-1 DO master[2000:2511,2000:2511,i]=mrdfits(master_files[i])

;calculate DEM
print,'Calculating DEM'
aia_firdem_wrapper,filters,exp_times,master,dem
print,'Saving DEM'
save,filename=strcompress('dem_'+files+'.sav',/remove_all),/all

;set temperature range for DEM temperature image
logtmin=6.1
logtmax=6.3

;plot DEM temperature image
print,'Plotting DEM images'
firdem_emwtemp_image_color,dem,logtmin,logtmax,image,emwmt=emwmt,emtot=emtot

;Save image output
print,'Saving DEM images'
save,image,filename=strcompress('DEMtemperature_'+files+'_'+string(logtmin)+$
                          '_'+string(logtmax)+'.sav',/remove_all)
;temperature
save,emwmt,filename=strcompress('emwmt_'+files+'.sav',/remove_all)
;emission measure
save,emtot,filename=strcompress('emtot_'+files+'.sav',/remove_all)

;plot DEM images
pos=[0.05,0.15,0.8,0.90]        ;Device coordinates for images

set_plot,'ps'
device,/encapsulated,/decomposed,/color,bits_per_pixel=8,$
  filename=strcompress('DEMtemperature_'+files+'_'+string(logtmin)+'_'+$
                       string(logtmax)+'.eps',/remove_all),$
  xsize=15,ysize=15,/cm
plot_image,image,true=3,xtickformat="(A1)",ytickformat="(A1)"
device,/close

device,/encapsulated,bits_per_pixel=8,$
  filename=strcompress('emwmt_'+files+'.eps',/remove_all),$
  xsize=15,ysize=15,/cm
plot_image,emwmt,xtickformat="(A1)",ytickformat="(A1)",position=pos
fsc_colorbar,range=[min(emwmt,/nan),max(emwmt,/nan)],$
  position=[pos(2),pos(1),pos(2)+0.05,pos(3)],format='(F4.2)',$
  title=textoIDL('log(T) [K]'),$
  /vertical,/right,charsize=0.9
device,/close

device,/encapsulated,bits_per_pixel=8,$
  filename=strcompress('emwmt_'+files+'_'+string(logtmin)+'_'+$
                       string(logtmax)+'.eps',/remove_all),$
  xsize=15,ysize=15,/cm
plot_image,emwmt,xtickformat="(A1)",ytickformat="(A1)",position=pos
fsc_colorbar,range=[logtmin,logtmax],$
  position=[pos(2),pos(1),pos(2)+0.05,pos(3)],format='(F4.2)',$
  title=textoIDL('log(T) [K]'),$
  /vertical,/right,charsize=0.9
device,/close

device,/encapsulated,bits_per_pixel=8,$
  filename=strcompress('emtot_'+files+'.eps',/remove_all),$
  xsize=15,ysize=15,/cm
plot_image,emtot,xtickformat="(A1)",ytickformat="(A1)",position=pos
fsc_colorbar,range=[min(emtot),max(emtot)],$
  position=[pos(2),pos(1),pos(2)+0.05,pos(3)],format='(E8.2)',$
  title=textoIDL('cm^{-5}'),$
  /vertical,/right,charsize=0.9
device,/close

set_plot,'x'

;;;;;Calculate Thermal Bremsstrahlung emission from DEM, White 2000;;;;;
;Define constants
flux_constant=9.78d-3           ;from White 2000, eqn 1
k_B=1.38d-16                    ;g cm^2 s^-2 K^-1
c=3d10                          ;cm s^-1
Nhe_Nh=1d-1                     ;He/H from White 2000
Nfe_Nh=1.56d-4                  ;Iron/H from White 2000
f=2.8d9                         ;2.8 GHz
AIA_pixel=double(0.59)          ;arcsec
arcsec_to_radian=double(2.*!pi/(360.*60.*60.)) ;radians per arcsecond
solid_angle=(AIA_pixel*arcsec_to_radian)^2 ;arcesc^2
cgs_to_sfu=(1d-22)^(-1)*(1d-7)*(1d-2)^(-2) ;[sfu]=[10^-22 W m^-2 Hz^-1]

constant_product=flux_constant*2*k_B*(c^(-2))*(1+4*Nhe_Nh)*solid_angle

;Create bremsstrahlung and T arrays
bremsstrahlung=dblarr(n_elements(dem.coffs[*,0,0]),$
                      n_elements(dem.coffs[0,*,0]))
T=1.*10^(dem.T0A)

;Create dT array (for discrete integration)
dT=fltarr(n_elements(T))
FOR i=0,n_elements(T)-1 DO BEGIN
    CASE i OF 
        0: dT[i]=T(i+1)-T(i)
        n_elements(T)-1: dT[i]=T(i)-T(i-1)
        ELSE: dt[i]=(0.5)*(T(i+1)-T(i-1))
    ENDCASE
        
ENDFOR

;Calculate integral of temperature from each pixel individually
print,'Calculating Bremsstrahlung flux'
FOR i=0,n_elements(dem.coffs[*,0,0])-1 DO BEGIN
    FOR j=0,n_elements(dem.coffs[0,*,0])-1 DO BEGIN
        bremsstrahlung[i,j]=total((T^(-0.5))*dem.coffs(i,j,*)*$
                                  (24.5-alog(T/f))*dT)
    ENDFOR
ENDFOR

;Scale temperature integral by the solid angle and flux constant
bremsstrahlung=bremsstrahlung*constant_product

;Save and plot predicted 10.7 cm radio flux from bremsstrahlung emission.
print,'Saving Bremsstrahlung image'
save,bremsstrahlung,filename=strcompress('bremsstrahlung_'+files+$
                                         '.sav',/remove_all)

print,'Total bremsstrahlung flux [erg cm^-2 s^-1 Hz^-1]:  '+$
  string(total(bremsstrahlung))
print,'sfu:  '+string(total(bremsstrahlung)*cgs_to_sfu)

set_plot,'ps'
device,/encapsulated,bits_per_pixel=8,$
  filename=strcompress('bremsstrahlung_'+files+$
                       '.eps',/remove_all),$
  xsize=15,ysize=15,/cm
plot_image,bremsstrahlung,xtickformat="(A1)",ytickformat="(A1)",position=pos
fsc_colorbar,range=[min(bremsstrahlung),max(bremsstrahlung)],$
  position=[pos(2),pos(1),pos(2)+0.05,pos(3)],format='(E8.2)',$
  title=textoIDL('erg cm^{-2} s^{-1} Hz^{-1}'),$
  /vertical,/right,charsize=0.9
device,/close
set_plot,'x'

stop

end