from pylab import *

def planck(lam, T):
    
    c= 2.99792458e8
    k= 1.380662e-23
    h= 6.626180e-34
    return 2*pi*h*c*c/(lam**5 * (exp(h*c/(lam*k*T))-1.))


# solar spectrum
solar_spectrum=loadtxt('/home/claudia/libRadtran/data/solar_flux/kurudz_1.0nm.dat')

d= 1.49e11      # sun-Earth distance
r= 1.3914e9/2.  # radius of the sun


# irradiance from 300 to 500 nm
sw_spectrum=loadtxt('irradiance.out')

figure(1, figsize=(10,7))

plot(sw_spectrum[:,0], sw_spectrum[:,1]+sw_spectrum[:,2], label='irradiance')
plot(solar_spectrum[:,0], solar_spectrum[:,1], 'k', label='TOA')
plot(sw_spectrum[:,0], planck(sw_spectrum[:,0]*1e-9, 6000)*1e-6*(r/d)**2, label='Planck')
xlim(200,2500)
title('irradiance spectrum at the surface')
xlabel('wavelength [nm]')
ylabel('irradiance [mW/ (m^2 nm)]')
legend(loc=1)
savefig('solar_irrad.png')

# transmittance from 300 to 500 nm
trans_spectrum=loadtxt('transmittance.out')

figure(2, figsize=(10,7))
plot(trans_spectrum[:,0], trans_spectrum[:,1]+trans_spectrum[:,2])
xlim(200,2500)
title('transmittance spectrum at the surface')
xlabel('wavelength [nm]')
ylabel('transmittance')
savefig('transmittance.png')