PRO readcf_id,line,init_model,id,cf,tau,h,prof,dlamax
; 
; Reading contribution functions from the database
; 
; INPUTS:
;   line - line name (string):
;       'La','Lb','Lc','Ld', 'Ha','Hb','Hc', 'Pa','Pb','Ba' for hydrogen
;       'CaK','CaH', '8498','8542','8662' for calcium II
;       'Mgk','Mgh', '2790.8','2797.9','2798.0' for magnesium II
;   init_model (string) in the form 'C' defines the folder name
;   id (integer) - ID number of the model in the folder 'Models/VAL_'+init_model
;   dlamax - maximum delta lambda [A] for plotting (optional)
; 
; OUTPUTS:
;   cf - contribution function, 2D FLOAT array [delta_lambda, height]
;   tau - optical depth in the line, 2D FLOAT array [delta_lambda, height]
;   h  - geometrical height in the model [km]
;   prof - synthetic profile of the line, DOUBLE(2,*) [delta_lambda, intensity]
;          units [A, erg/cm^2/s/sr/Hz], delta_lambda starts at 0 A
;
; CALLING: READ_MODEL_ID.PRO, READ_PROFILE_ID.PRO, CGLOADCT.PRO
; 
; EFFECT: Makes a contour plot from max(cf)/1000 to max(cf) in logarithmic scale,
;         together with the line profile (red) and tau=1 contour (blue).
; 
on_error,1

id=id<2988>0

; LINES DATA
name=['La','Lb','Lc','Ld','Ha','Hb','Hc','Pa','Pb','Ba', $
      'CaK','CaH','8662','8498','8542','Mgh','Mgk','2790.8','2798.0','2797.9']
elem=['HYD','HYD','HYD','HYD','HYD','HYD','HYD','HYD','HYD','HYD', $
      'CA','CA','CA','CA','CA','MG','MG','MG','MG','MG']
ll=  [15,36,57,77,97,137,165,193,221,249, $
      8,36,64,82,100,26,66,106,138,170]
ul=  [35,56,76,96,136,164,192,220,248,276, $
      35,63,81,99,123,65,105,137,169,201]

w=where(line eq name)
if w eq -1 then message,'Unexpected line name encountered'
w=w[0]

filename ='Contrib_Functions/'+'cf_VAL'+init_model+elem[w]+'/conf'+string(id,format='(I05)')+'.DAT'

READ_MODEL_ID,init_model,id,modin,flag
if flag eq 0 then message,'Contribution function not found (NaN)'
h=reform(modin[0,1:102])
READ_PROFILE_ID,line,init_model,id,prof

nlam=ul[w]-ll[w]+1     ; number of wavelengths for line name[w]
cf=fltarr(102,nlam)    ; contribution function array
tau=fltarr(102,nlam)   ; optical depth array

openr,1,filename
; Skip everything before the line name[w]
for l=0,102L*ll[w]-1 do readf,1,tau1,cf1
; Read the data of line name[w]
for j=0,nlam-1 do begin    ; loop over lambda
  for i=0,101 do begin     ; loop over depth
    readf,1,tau1,cf1
    tau[i,j]=tau1
    cf[i,j]=cf1
  endfor
endfor
close,1

cf=transpose(cf)    ; horizontal: lambda, vertical: depth
tau=transpose(tau)

; Plotting
window,1,xs=800,ys=600
CGLOADCT,3,/reverse

if n_params() lt 8 then begin
  sip=size(prof)
  dlamax=prof[0,sip[2]-1]
endif
xran=[0,dlamax]       ; exact wavelength range

; Contour levels
; levels=[1.e-15,3.3e-15,6.7e-15,1.e-14,3.3e-14,6.7e-14,1.e-13]

cfmin=alog10(max(cf)/1000.)          ; logarithmic contour levels
cfmax=alog10(max(cf))                ; from max(cf)/1000 to max(cf)
levstep=(cfmax-cfmin)/10.
levels=cfmin+findgen(11)*levstep

;contour,cf,prof[0,*],h,lev=levels,xra=xran,yra=[0,2500],/xst,/yst, $
contour,alog10(cf>cfmin),prof[0,*],h,lev=levels,xra=xran,yra=[-100,2500],/xst,/yst, $
  title='Line '+line+'  model '+init_model+'-'+string(id,format='(I04)'), $
  xtit='!4Dk!3 ['+string("305B)+']',ytit='Height [km]',charsiz=1.5,/fill

loadct,1,/silent  
contour,tau,prof[0,*],h,lev=[0,1],xra=xran,yra=[-100,2500],xst=5,yst=5,/noerase, $
  col=150,charsiz=1.5                ; blue line

loadct,3,/silent
plot,prof[0,*],prof[1,*],xra=xran,xst=5,yst=5,/noerase,col=150,charsiz=1.5
                                     ; red line

END