# from Clemens' flux calc &  Robin's paper [Owen et al, JSR 16, 2009, Silicon PIN diodes]
#--------- formula --------------------------------------------
# flux = I_ph * eps_Si /e/E/( 1-exp(A_Si*rsi*t_Si ) ) 
#     with Al and air atten:
# flux= flux * exp(A_Al*ral*t_Al) * exp(A_air*rair*t_air)
# 
# with A_ being the photoelectric cross sections of the materials 

#--------- parameters -----------------------------------------
# density of Si
rsi = 2.33      # g/cm^3
# density of Al 
ral = 2.699     # g/cm^3
# density of air 
rair    =       1.205e-3  # g/cm^3

eps_si = 3.62 # eV energy req. for charge separation in Si (generation of el/hole pairs) 
e      = 1.602e-19    # As , elementary charge 

t_si = 0.0012 # thickness of diode in cm  (== 12 micron) ; [* 10000] 
t_al = 0.002 # thickness of diode in cm ; (== 20 micron)

#--------- input---------------------------------------------------

cur = float(get_string('Please enter the measured diode current [in mA]', 0.1))
ep = float(get_string('Please enter the photon energy  [in keV]', 12.4 ))
t_si= float(get_string('Please enter the thickness of the Si layer  [in micron]', t_si))
t_si = t_si/10000.  # --> cm 
t_al = float(get_string('Please enter the thickness of the Al layer  [in micron, 0 if not available]', t_al))
t_al = t_al/10000
t_air = float(get_string('Please enter the pathway in air  [in mm]', 165.+100.))


#--------- calc----------------------------------------------------
# energy deposit in Silicon

polys=[4.158,- 2.238, - 0.477, 0.0789]

hlp_si= poly(math.log10(ep), polys)
A_Si= 10.0 ** hlp_si
efact=math.exp(-A_Si*rsi*t_si)

sifact=1.0-efact

#help, cur, eps_si, ep

fl0=cur* eps_si/1.602/ep/sifact *1.e13 


# Aluminium attenuation

# ratio of photoelectric cross section to density for Aluminium
polyal= [4.106, - 2.349, - 0.413, 0.0638 ]
hlp_al= poly(math.log10(ep), polyal)
A_Al= 10.0 ** hlp_al

# attenuation due to aluminium
alfact=math.exp(-A_Al*ral*t_al)

# Air attenuation

# ratio of photoelectric cross section to density for air
polyair= [3.153, - 1.026, - 2.348,  0.928] 
hlp_air= poly(math.log10(ep), polyair)
A_air= 10.0 ** hlp_air

# attenuation due to air 
airfact=math.exp(-A_air*rair*t_air/10.)

# total flux from photocurrent 

fl=fl0/alfact/airfact

f = fl 

msg =      '  Energy:       ' + '%7.4f' % e + ' keV\n'
msg = msg+ '  Diodecurrent: ' + str(cur).strip() +' mA\n'
msg = msg+ '\n'
msg = msg+ '  Thickness of active Si layer:            ' + str(t_si*10000.).strip() + ' micron\n'
msg = msg+ '  Thickness of Al layer in front of diode: ' + str(t_al*10000.).strip() + ' micron\n'
msg = msg+ '  Length of path in air in front of diode: ' + str(t_air).strip() + ' mm\n'
msg = msg+ '\n'
msg = msg+ '  ===> flux:         %8.2E  photons / s' % f

print msg
show_message(msg, "flux_diode", False)