pro opacity_read,wave,opacity

; v1 - 6jan09 - dhc - reads in cmfgen kappa model from Janos
;                     plots opacity vs. lamba for specified radius

; the above comment refers to opacity_plot.pro (v1) on which 
; this routine is based

; v1 - 7jan09 - dhc - reads in cmfgen opacity model and passes wavelength 
;                     and opacity arrays


rstar=123.6                ;Rstar in units of 1e10 cm

; READING INPUT FILE


;openr,unit,'opacity.dat',/get_lun

; OPACITY MODEL SENT BY JANOS DEC 2008
;openr,unit,'opacities/CONT_OPAC.Mdot=3E-6',/get_lun

; OPACITY MODEL SENT BY JANOS FEB 2009 - incl. Mg and Si
openr,unit,'opacities/CONT_OPAC.Mdot=1.7E-6.feb09',/get_lun
print, 'USING OPACITY MODEL: opacities/CONT_OPAC.Mdot=1.7E-6.feb09'

header=''
num_rad=5L
num_wave=5L

;read the first line, pull out the number of depth points
readf,unit,num_rad,header
print,num_rad
print,header

;read the second line, pull out the number of frequency points
readf,unit,num_wave,header
print,num_wave
print,header

; RADIUS GRID

;read the blank line
readf,unit,header
print,header

;read the title block of the radius grid
readf,unit,header
print,header

;read the 6 columns by 10 rows block of radii -- this 6x10 will change if num_rad changes
;but this code won't automatically account for that now

radii=fltarr(60)
data=fltarr(6,10)
readf,unit,data
print,data

for i = 0,9 do begin
for j = 0,5 do begin

radii[i*6+j] = data[5-j,9-i]

endfor
endfor

print,' '
print,'radii array 0 through 59'
print,radii/rstar

; VELOCITY GRID - USING SAME TECHNIQUE AS FOR RADIUS GRID

;read the blank line
readf,unit,header
print,header

;read the title block of the velocity grid
readf,unit,header
print,header

;read the 6 columns by 10 rows block of velocities -- this 6x10 will change if num_rad changes
;but this code won't automatically account for that now

velocities=fltarr(60)
data=fltarr(6,10)
readf,unit,data
print,data

for i = 0,9 do begin
for j = 0,5 do begin

velocities[i*6+j] = data[5-j,9-i]

endfor
endfor

print,' '
print,'velocities array 0 through 59'
print,velocities


; DENSITY GRID - USING SAME TECHNIQUE AS FOR RADIUS GRID

;read the blank line
readf,unit,header
print,header

;read the title block of the density grid
readf,unit,header
print,header

;read the 6 columns by 10 rows block of densities -- this 6x10 will change if num_rad changes
;but this code won't automatically account for that now

densities=fltarr(60)
data=fltarr(6,10)
readf,unit,data
print,data

for i = 0,9 do begin
for j = 0,5 do begin

densities[i*6+j] = data[5-j,9-i]

endfor
endfor

print,' '
print,'densities array 0 through 59'
print,densities


; CLUMPING FACTOR GRID - USING SAME TECHNIQUE AS FOR RADIUS GRID

;read the blank line
readf,unit,header
print,header

;read the title block of the clumping factor grid
readf,unit,header
print,header

;read the 6 columns by 10 rows block of clumping factors -- this 6x10 will change if num_rad changes
;but this code won't automatically account for that now

clumpfac=fltarr(60)
data=fltarr(6,10)
readf,unit,data
print,data

for i = 0,9 do begin
for j = 0,5 do begin

clumpfac[i*6+j] = data[5-j,9-i]

endfor
endfor

print,' '
print,'clumping factors array 0 through 59'
print,clumpfac

; PLOTTING CLUMPING FACTOR VS. RADIUS

;window
;!p.font=0
;plot,radii/rstar,clumpfac,xrange=[1.,10.]

; NOTE: F_C = 0.15 PRETTY MUCH EVERYWHERE, EXCEPT RIGHT AT THE BASE


; READING IN THE 10,000 FREQUENCY BLOCKS OF 60 OPACITY VALUES

; first, handle the blank line and header of the kappa table

;read the blank line
readf,unit,header
print,header

;read the title block of the entire kappa table
readf,unit,header
print,header



; NOW, BELOW HERE, IS THE FIRST BLOCK OF ACTUAL OPACITY DATA
; ok, CHECKS OUT
; NEXT, LOOP OVER THE FIRST THREE BLOCKS 
; ok, CHECKS OUT
; NEXT, LOOP OVER ALL NUM_WAVE BLOCKS

print,'num_wave= ',num_wave
wave=fltarr(num_wave)
opacities=fltarr(num_wave,num_rad)

; ACTUALLY, FIRST PRINT OUT THE NUM_RAD RADII IN TERMS OF RSTAR 
; IN A COLUMN

for kk = 0,num_rad-1 do begin
print,kk,radii[kk]/rstar
endfor

; NOW, THE USER CAN CHOOSE WHICH RADIUS TO PICK

plot_index=2L
print,'which radial index do you want to choose the opacity for?'
read,plot_index
print,plot_index

for k = 0,num_wave-1 do begin


;read the blank line
readf,unit,header

;read the wavelength of the ith block of the kappa table
readf,unit,wavetmp
wave[k]=wavetmp

;read the 6 columns by 10 rows block of densities -- this 6x10 will change if num_rad changes
;but this code won't automatically account for that now

data=fltarr(6,10)
readf,unit,data

for i = 0,9 do begin
for j = 0,5 do begin

opacities[k,i*6+j] = data[5-j,9-i]

endfor
endfor


endfor


; CLOSING FILE

free_lun,unit


; NOW, PRINT OUT THE USER-SPECIFIED INDEX AND RADIUS, VELOCITY, DENSITY, F_C

print,'you have specified a plot for index: ',plot_index
print,'at radius = ',radii[plot_index]/rstar
print,'where the values of velocity, density, clumping factor are:'
print,velocities[plot_index],densities[plot_index],clumpfac[plot_index]
print,' '

; AND SPOT CHECK A FEW OF THE OPACITIES

;print,'the opacities at this radius, for a few specific wavelengths are:'
;print,'wave[0],opacities[0,plot_index]'
;print,wave[0],opacities[0,plot_index]
;print,'wave[9999],opacities[9999,plot_index]'
;print,wave[9999],opacities[9999,plot_index]

; AND PUT THE USER-SPECIFIED RADIAL SLICE OF THE 2-D OPACITIES ARRAY INTO A 1-D ARRAY

opacity=fltarr(num_wave)
opacity=opacities(*,plot_index)
;print,wave[0],opacity[0]
;print,wave[9999],opacity[9999]


return
end