;Photodissociation and photoionization graphs
;Employer: Ewine van Dishoek, Sterrewacht Leiden
;Author: Sander de Kievit
;Use with care, blatant errors will and might still be present
;Needs adjustments to work on other systems than strw.leidenuniv.nl
;WARNING: READS, DELETES AND RECREATES FILES ON YOUR FILE SYSTEM!!!
;WARNING: ALPHA RELEASE

PRO photo
	FORWARD_FUNCTION splits_arrays
	COMMON arrays, angstrom, cross_section, molecules, filenames, Tbb, dilfac, angstrom_l, cross_l
	COMMON strings, filename, photo, title, mol, graphicsname, molecule 
	COMMON floats, k_line_draine, k_cont_draine, k_line_bb, k_cont_bb
	COMMON constants, kb, c, h, pi
	COMMON luns, drainelun, bblun, mollun
	COMMON formats, format
	COMMON store, k_line_bb_nh3tonh, k_line_draine_nh3tonh, k_line_bb_nh3, k_line_draine_nh3, k_cont_bb_h2otoo, k_cont_draine_h2otoo, k_cont_bb_h2o, k_cont_draine_h2o

	format = '(E8.1)'

;	define the physical constants
	kb = 1.380658D-16
	h = 6.626176D-27
	c = 2.99792458D10
	pi = 3.141593D0

;	Indicate which blackbody temperatures we will be using (array - vector)
;	and initialize the dilution factor array (vector)
;	Tbb=double([4000.0, 5000.0, 10000.0, 20000.0, 50000.0])
	Tbb=double([4000.0, 6000.0, 10000.0, 20000.0])
	dilfac=DBLARR(size(Tbb, /N_ELEMENTS))

	;Load colour table, for future use maybe?
	loadct,13

	;Create string array with filenames (pd = dissociation, pi = ionization)
	p_ion=findfile('pd/*.pi')
	p_dis=findfile('pd/*.pd')

	;Read string array from file for converting the file name to the name of the molecule (2 x n array)
	readcol, "molecules", filenames, molecules, FORMAT='A,A', DELIMITER=',' ,/SILENT

	;Determine number of objects in the arrays
	p_ion_aantal=size(p_ion,/N_ELEMENTS)
	p_dis_aantal=size(p_dis,/N_ELEMENTS)
;	p_ion_aantal = 1
;	p_dis_aantal = 1
	;Set the plot ('ps') is for file, 'X' is for xserver and set the font used in the plots.
	set_plot, 'ps'
	!P.FONT = 0 

	;Open the file where we will put the dissociation rates for the draine field and the temperatures fields
	OPENW, drainelun, "/disks/strw7/dekievit/werk/rates-draine", /GET_LUN
	OPENW, bblun, "/disks/strw7/dekievit/werk/rates-bb", /GET_LUN

	;Initialize the blackbody rates file with the numer of different temperatures
	printf, bblun, strcompress(string(size(Tbb, /N_ELEMENTS)), /REMOVE_ALL)

	;Give the temperatures with the first element in the row indicating these are the temperatures (Temp)
	Tbbstring = "Temp"
	FOR i=0, size(Tbb, /N_ELEMENTS) - 1 DO Tbbstring = Tbbstring + "," + string(Tbb[i])
	printf, bblun, strcompress(Tbbstring, /REMOVE_ALL)

	; Calculate the dilution factors corresponding to the bb temperatures
	dilfac = calc_dil(Tbb)
	print, dilfac
	
	;Set the string to photoionization or photodissocation
	;This string appears at the top of the plots.
	photo = "Photoionization "
	make_plots, p_ion_aantal, p_ion
	photo = "Photodissociation "
	make_plots, p_dis_aantal, p_dis
	FREE_LUN, drainelun
	FREE_LUN, bblun
END


;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;----------------------------- FUNCTIONS --------------------------------------
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------

FUNCTION calc_dil, Tbb
	COMMON constants
	device, filename='draine-bb.ps', /color
	device, /TIMES
	device, /ISOLATIN1
;	Initialize an array with wavelengths in Angstrom (from 912 A to 2050 A)
;	To be consistent with Faraday article, use step sizes of 20 l.
;	wavelength = DINDGEN(58)
;	wavelength = (912.00 + wavelength*20.)
	wavelength = DINDGEN(1139)
	wavelength = (912.00 + wavelength)
;	Calculate the draine field
	I_draine = double(h*c/(wavelength*1.0D-8))*calc_draine(wavelength)
;	Initialize the dilution array
	dilution_array = DBLARR(size(Tbb,/N_ELEMENTS))
	FOR num=0, size(Tbb, /N_ELEMENTS)-1 DO BEGIN
		T = double(Tbb[num])
;	Implementation of draine field in the reference on: "http://adsabs.harvard.edu/abs/1988rcac.book...49V"
;	Implement blackbody formula from Rybicki and Lightman p. 22 eq. 1.52 integrated over the entire sphere (4*pi)
		Bb = double((4.0*pi*2.0*h*c^(2.0))/((wavelength*1.D-8)^(5.0)))/double(exp(double((h*c)/(wavelength*1.D-8*kb*T))) - 1.0)
;	Calculate the dilution factors for normalization of the black body to the draine field.
		dilution_array[num] = double(tsum(wavelength, I_draine)) / double(tsum(wavelength*1.D-8, Bb))
		print, T, dilution_array[num]
	ENDFOR

;	Do the samething, but then for the figure (note that the units and scales can be different!!)
	wavelength = DINDGEN(2089)
	wavelength = (912.00 + wavelength)
	I_draine = calc_draine(wavelength)
	plot, wavelength, I_draine, title="Draine field vs scaled blackbody fields", $
			xtitle="Wavelength (" + STRING(197B) + ")", ytitle="Flux (photon cm!E-2!N s!E-1!N " + STRING(197B) + "!E-1!N)", /YLOG, $
			yrange=[1.D2,1.D8], xrange=[800,3000], xminor=1, xstyle=1 
	newcolor = INTARR(size(Tbb,/N_ELEMENTS))
	FOR num=0, size(Tbb, /N_ELEMENTS)-1 DO BEGIN
		T = double(Tbb[num])
;	Implementation of draine field in the reference on: "http://adsabs.harvard.edu/abs/1988rcac.book...49V"
;	Implement blackbody formula from Rybicki and Lightman p. 22 eq. 1.52 integrated over the entire sphere (4*pi)
		Bb = double((4.0*pi*2.0*c)/((wavelength*1.D-8)^(4.0)))/double(exp(double((h*c)/(wavelength*1.D-8*kb*T))) - 1.0)
		newcolor[num] = round((T/max(Tbb))*256.0)
		print, min(Bb*dilution_array[num]*1.D-8), max(Bb*dilution_array[num]*1.D-8)
		oplot, wavelength, Bb*dilution_array[num]*1.D-8, color=newcolor[num]
	ENDFOR

	legend, ['T=' + STRING(strcompress(round(Tbb),/REMOVE_ALL))], line=0, color=newcolor
	device, /close
	RETURN, dilution_array
END

PRO make_plots, array_aantal, array_filenames
	COMMON arrays
	COMMON strings
	COMMON floats
	COMMON luns
	COMMON formats
	COMMON store

	FOR nummer=0,(array_aantal-1) DO BEGIN
		; initialize string for in the plots
		; obtain the filename
		filename=array_filenames[nummer]
		; if we don't want this photodissociation, say so here!
;		if (filename EQ "pd/h2o.pd" || filename EQ "pd/nh3.pd") THEN print, "Found one!"
;		ENDIF ELSE BEGIN
		; Open for read, read the first line into two strings
		OPENR, lun, filename, /GET_LUN
		header = STRARR(2)
		READF, lun, header
		rows_line = 0  ; Number of rows in the line spectrum (initialize to 0)
		; Read the correct information from the file
		READS, header(1), rows_line, FORMAT='(A17, x, I0)'
		; Move to the position right after the line data ends:
		header = STRARR(rows_line+1)
		; read all lines of the line spectrum and the line indicating how many
		; lines of data we have of continuum data:
		READF, lun, header
		rows_continuum = 0 ; Number of rows in the continuum spectrum
		; Determine the number of lines that the continuum data has
		; Don't use the other data read (will be done later with readcol)
		READS, header(rows_line), rows_continuum, FORMAT='(A17, x, I0)'
		; Free the logical file unit
		FREE_LUN, lun
		; obtain the molecule or atom from the filename
		mol = STRMID(filename, 3, STRPOS(filename, ".p") - 3)

		; IDL lacks normal regular expressions, so a bit of nasty code
		; to convert things to upper and lowerscript
		molecule = molecules[WHERE(filenames EQ mol)]
		molecule = REPSTR(molecule, '2', '!I2!N')
		molecule = REPSTR(molecule, '3', '!I3!N')
		molecule = REPSTR(molecule, '4', '!I4!N')
		molecule = REPSTR(molecule, '6', '!I6!N')
		molecule = REPSTR(molecule, '+', '!E+!N')
		molecule = REPSTR(molecule, '-', '!E-!N')
		molecule = REPSTR(molecule, '=', '-')
		if (filename EQ "pd/h2o.pd") THEN molecule = molecule + " to OH"
		if (filename EQ "pd/nh3.pd") THEN molecule = molecule + " to NH!I2!N"

		; Inlezen van kolommen van de verschillende datafiles
		readcol, filename, number_array, angstrom, cross_section, /SILENT

		; Initialize to 0 to have easier testing conditions later!
		k_line_draine = 0.0
		k_cont_draine = 0.0

		; Create the arrays where the rates (k) can be read into
		k_line_bb = DBLARR(size(Tbb, /N_ELEMENTS))
		k_cont_bb = DBLARR(size(Tbb, /N_ELEMENTS))

		; Create the filename for the graphics file (molecule_pd/pi.ps)
		graphicsname = mol + "_" + STRMID(filename, STRPOS(filename,".p")+1, STRPOS(filename,".p") + 2)
	
		; Create the title of the plots of line data
		title = photo + "line channels of "
		hc3h=0
		if (filename EQ "pd/hc3h.pd") THEN hc3h=1
		if (hc3h) THEN print, "HC3H!!", rows_line, rows_continuum

		; If we have line data, make a stik_diagram
		IF rows_line NE 0 THEN BEGIN
			stik_diagram, rows_line
			k_line_bb = bb(rows_line, angstrom_l, cross_l, Tbb, dilfac, '-line')
			if (hc3h) THEN print, "HC3H!! stikdiagram"
		ENDIF
		; Create the title of the plots of continuum data
		title = photo + "continuum data of "

		; If we have continuum data, make the plot
		IF rows_continuum NE 0 THEN BEGIN
			spline_interpol_diagram, rows_line, rows_continuum
			k_cont_bb = bb(size(angstrom, /N_ELEMENTS), angstrom, cross_section, Tbb, dilfac, '-cont')
		ENDIF
		; If we have both line and continuum data, also include a total in the file
		IF (k_line_draine NE 0.0 && k_cont_draine NE 0.0) THEN BEGIN
			data = graphicsname + "-total," + string(k_line_draine + k_cont_draine, FORMAT=format)
			printf, drainelun, strcompress(data, /REMOVE_ALL)
		ENDIF
		IF (ARRAY_EQUAL(k_line_bb, 0.0) || ARRAY_EQUAL(k_cont_bb, 0.0)) THEN BEGIN
		ENDIF ELSE BEGIN
			kstring = graphicsname + "-total"
			FOR num=0, size(Tbb, /N_ELEMENTS)-1 DO BEGIN
				kstring = kstring + "," + string(k_line_bb[num]+k_cont_bb[num], FORMAT=format)
			ENDFOR
			printf, bblun, strcompress(kstring, /REMOVE_ALL)
		ENDELSE
		
;		if filename EQ "pd/nh3tonh.pd" THEN BEGIN
;		ENDIF if filename EQ "pd/h2otoo.pd" THEN BEGIN
		
		if filename EQ "pd/nh3tonh.pd" THEN BEGIN
			k_line_bb_nh3tonh=k_line_bb
			k_line_draine_nh3tonh=k_line_draine
		ENDIF
		if filename EQ "pd/nh3.pd" THEN BEGIN
			k_line_bb_nh3=k_line_bb
			k_line_draine_nh3=k_line_draine
		ENDIF
		if filename EQ "pd/h2otoo.pd" THEN BEGIN
			k_cont_bb_h2otoo=k_cont_bb
			k_cont_draine_h2otoo=k_cont_draine
		ENDIF
		if filename EQ "pd/h2o.pd" THEN BEGIN
			k_cont_bb_h2o=k_cont_bb
			k_cont_draine_h2o=k_cont_draine
		ENDIF	
		OPENW, mollun, graphicsname + ".dat", /GET_LUN
		printf, mollun, mol + "  " + photo
		printf, mollun, "Draine " , Tbb
		printf, mollun, k_line_draine, k_line_bb
		printf, mollun, k_cont_draine, k_cont_bb
		FREE_LUN, mollun
	ENDFOR
	IF photo EQ "Photodissociation " THEN BEGIN
		kstring_nh3 = "nh3_pd-all"
		kstring_h2o = "h2o_pd-all"
		data = "nh3_pd-all," + string(k_line_draine_nh3tonh + k_line_draine_nh3, FORMAT=format)
		printf, drainelun, strcompress(data, /REMOVE_ALL)
		data = "h2o_pd-all, " + string(k_cont_draine_h2otoo + k_cont_draine_h2o, FORMAT=format)
		printf, drainelun, strcompress(data, /REMOVE_ALL)
		FOR num=0, size(Tbb, /N_ELEMENTS)-1 DO BEGIN
			kstring_nh3 = kstring_nh3 + "," + string(k_line_bb_nh3tonh[num]+k_line_bb_nh3[num], FORMAT=format)
			kstring_h2o = kstring_h2o + "," + string(k_cont_bb_h2otoo[num]+k_cont_bb_h2o[num], FORMAT=format)
		ENDFOR
		printf, bblun, strcompress(kstring_nh3,/REMOVE_ALL)
		printf, bblun, strcompress(kstring_h2o,/REMOVE_ALL)
	ENDIF
END


PRO stik_diagram, rows
	COMMON arrays
	COMMON strings
	COMMON floats
	COMMON formats
	angstrom_l = angstrom[0:rows-1]
	cross_l = cross_section[0:rows-1]
	device, filename=graphicsname + '-line.ps', /color
	device, /TIMES
	IF (MAX(angstrom_l)-MIN(angstrom_l) LT 400) $
		THEN xrange = [MIN(angstrom_l) - 100, MAX(angstrom_l) + 100] $
	ELSE xrange = [MIN(angstrom_l) - 50, MAX(angstrom_l) + 50]
		PLOT, angstrom_l, cross_section[0:rows-1]*10.0^(15), PSYM=3, SYMSIZE=0.1, $
		   title=title + molecule, xtitle="Wavelength (" + STRING(197B) + ")", ytitle="Cross section (x 10!E-15!N cm!E2!N " + STRING(197B) +")", $
		   CHARSIZE=1.3, CHARTHICK=4, XMARGIN=[12,3], XTHICK=4, YTHICK=4, XRANGE = xrange
	FOR i=0, rows-1 DO BEGIN
		PLOTS, angstrom[i], 0
		PLOTS, angstrom[i], cross_section[i]*10.0^(15), /CONTINUE, THICK=4
	ENDFOR
	; Calculate photodissociation rate:
	k_line_draine = draine(rows, angstrom_l, cross_l, '-line')
	device, /close
END

PRO spline_interpol_diagram, rows, rows_continuum
	total_rows = rows + rows_continuum
	COMMON arrays
	COMMON strings
	COMMON floats
	COMMON constants
	COMMON luns
	cross_section = cross_section[rows:total_rows-1]
	angstrom = angstrom[rows:total_rows-1]
	total_rows = rows_continuum
	angstrom_lim=float([0.0, 0.0])
  ; Cut off the final part of the arrays if that is zero (Spline interpolation will fail miserably, if left like they are)
	WHILE cross_section[total_rows-1] EQ 0 DO BEGIN
		cross_section = cross_section[0:total_rows-2]
		angstrom_lim[1] = angstrom[total_rows-1]
		angstrom = angstrom[0:total_rows-2]
		total_rows = total_rows - 1
	ENDWHILE
	IF cross_section[0] EQ 0 THEN BEGIN
		cross_section = cross_section[1:total_rows-1]
		angstrom_lim[0] = angstrom[0];
		angstrom = angstrom[1:total_rows-1]
	ENDIF
	IF min(where(cross_section EQ 0)) NE -1 THEN BEGIN
		print, graphicsname 
	ENDIF
	IF (size(angstrom,/N_ELEMENTS) GT 1) THEN BEGIN
		device, filename=graphicsname + '-cont.ps', /color
		device, /TIMES
		SPLINE_P, double(angstrom), double(cross_section), XR, YR, tan0=[1,0]
		IF (min(YR) LT 0.0) THEN YR[WHERE(YR LT 0.0)] = 0.0
		plot, XR, YR*10.0^(18), title= title + molecule, xtitle="Wavelength (" + STRING(197B) + ")", $
		   ytitle="Cross section (x 10!E-18!N cm!E2!N)", $
		   CHARSIZE=1.3, CHARTHICK=4, XMARGIN=[7,3], XTHICK=4, YTHICK=4, THICK=4, YSTYLE=2
		oplot, angstrom, cross_section*10.0^(18), psym=4, THICK=4
		oplot, angstrom_lim, float([0.0, 0.0]), psym=4, color=256, THICK=4
		;if (angstrom_lim[0] ne 0.) then print, angstrom_lim, molecule
		k_cont_draine = draine(size(angstrom,/N_ELEMENTS), double(angstrom), double(cross_section), '-cont')
		device, /close
	ENDIF
END

FUNCTION calc_draine, angstrom
	I = DBLARR(size(angstrom, /N_ELEMENTS))
	I = (3.202833D15)/(angstrom^(3.0)) - (5.154206D18)/(angstrom^(4.0)) + (2.054625D21)/(angstrom^(5.0))
	IF (max(angstrom) GT 2000.00) THEN $ 
		I[WHERE(angstrom GT 2000.00)] = (7.315D2)*angstrom[WHERE(angstrom GT 2000.00)]^(0.7)
	IF (min(angstrom) LT 911.75) THEN $ 
		I[WHERE(angstrom LT 911.75)] = 0
	RETURN, I
END 

FUNCTION draine, rows, angstrom, cross, fileadd
	COMMON luns
	COMMON strings
	COMMON formats
	I = double(calc_draine(angstrom)*cross)
	IF (fileadd EQ "-line") THEN BEGIN
	;	draw_stick_diagram, angstrom, I, title + molecule + " in draine field", graphicsname + "-draine"
		k=TOTAL(I)
	ENDIF
	IF (fileadd EQ "-cont") THEN BEGIN
	;	draw_cont_diagram, angstrom, I, title + molecule + " in draine field", graphicsname + "-draine"
		k = tsum(angstrom, I) 
	ENDIF
	data = graphicsname + fileadd + "," + string(k, FORMAT=format)
	printf, drainelun, strcompress(data, /REMOVE_ALL)
	RETURN, k
END

FUNCTION bb, rows, angstrom, cross_section, Tbb, dilfac, fileadd
	COMMON strings
	COMMON constants
	COMMON luns
	COMMON formats
;	print, angstrom
;	print, cross_section
	kstring = ""
	k = DBLARR(size(Tbb, /N_ELEMENTS))
	FOR num=0, size(Tbb, /N_ELEMENTS)-1 DO BEGIN
		T = double(Tbb[num])
	; Create a string with the temperature for later use.
		tstring = strcompress(round(T),/REMOVE_ALL)
	; Create the array for the intensity (bb)
	; Convert to doubles (to avoid floating underflow exceptions)
		I = DBLARR(rows)
		wavelength = double(angstrom*1.D-8)
	; Convert to double for consistancy
		cross = double(cross_section)
	; Use formula 1.52 from Rybicki and Lightman (B_l(T)=(2hc^2/l^5)/(exp(hc/lkT)-1) (erg /(s cm^3 sr) )
;		Bb = double((4.0*pi*2.0*h*c^(2.0))/((wavelength*1.D-8)^(5.0)))/double((2.71828^((h*c)/(wavelength*1.D-8*kb*T)) - 1.0))
		I = (double((4.0*pi*2.0*c)/((wavelength)^(4.0))))/double(exp(double((h*c)/(wavelength*kb*T))) - 1.0)
	; Go over the normalization
		I = I * dilfac[num]
	; Multiply with the crosssection (Integrand)
	; Uncomment the below to check for the calculational mistakes (mostly of the order of E-16).
	;	IF (fileadd EQ "-cont") THEN print, norm, (tsum(angstrom, I_draine) - tsum(angstrom, I))/Int_draine $
	;	ELSE print, norm, (TOTAL(I_draine) - TOTAL(I))/Int_draine
		I = double(I * cross)
	; Do the integration, method depends on whether we are calculating continuous spectra or line emission
		IF (fileadd EQ "-cont") THEN BEGIN
;			draw_cont_diagram, angstrom, I, title + molecule + " in " + tstring + "K field", graphicsname + "-" + tstring
			k[num] = tsum(wavelength, I) 
		ENDIF ELSE BEGIN 
			k[num] = TOTAL(double(I*1.0D-8))
;			draw_stick_diagram, angstrom, I, title + molecule + " in " + tstring + "K field", graphicsname + "-" + tstring 
		ENDELSE
		kstring = kstring + "," + string(k[num], FORMAT=format)
	; Put the information into a file to be read out later by the script on the webserver
	ENDFOR
	data = graphicsname + fileadd + kstring
	printf, bblun, strcompress(data, /REMOVE_ALL)
	RETURN, k
END