PRO read_profile_id,line,init_model,id,prof,plot=plot
;
; Reading synthetic profiles 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
;
; OUTPUT:
;   prof - synthetic profile array DOUBLE (2,n) [delta_lambda, intensity]
;          units [A, erg/cm^2/s/sr/Hz], delta_lambda starts at 0 A
;
; KEYWORD: /plot - if set, a plot of the profile is made
;
;on_error,1

id=id<2988>0
tp=-0.03+(id/49)*0.001   ; temperature parameter
mp=0.2+(id mod 49)*0.1   ; column mass parameter
print,'Model T-par, m-par, ID:  ',tp,mp,id

; 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'

filename='Profiles/'+elem[w]+'/'+elem[w]+'_'+init_model+'/'+elem[w]+string(id,format='(I05)')+'.DAT'
n=ul[w[0]]
prof=fltarr(2,n+1)

openr,1,filename
for i=0,n do begin
  if elem[w] eq 'HYD' then readf,1,dum,dlam,intensity else readf,1,dlam,intensity
  prof[0,i]=dlam
  prof[1,i]=intensity
endfor
close,1

; Presence of NaNs
if total(finite(prof)) ne n_elements(prof) then message,'Profile not found (NaN)'

prof=prof[*,ll[w[0]]:ul[w[0]]]

;------------------ for MgII only ---------------------
;ma=max(prof[1,*],p_ma)
;print,'Position of maximum:',p_ma,prof[0,p_ma],' A'
;mi=min(prof[1,*],p_mi)
;mm=min(abs((ma-mi)/2.+mi - prof[1,p_ma:p_mi]),p)
;print,'Emission-peak HWHM:',p+p_ma,prof[0,p+p_ma],' A'
;------------------------------------------------------

if keyword_set(plot) then $
plot,prof[0,*],prof[1,*],xtit='!4Dk!3 ['+string("305B)+']', $
  ytit='I [erg cm!U-2!N s!U-1!N sr!U-1!N Hz!U-1!N]', $
  title='Line '+line+'  model '+init_model+'-'+string(id,format='(I04)'), $
  charsiz=1.5,xmargin=[14,4],col=0,back=255

;------------------ for MgII only -----------------
;plots,[prof[0,p+p_ma],prof[0,p+p_ma]],[0,ma],col=0
;--------------------------------------------------

END