closeall;  clear; 
run; 


################################################
# get simulation results

f_sim=getdata("T","f");      		
lambda_sim=c/f_sim;
Pabs = transmission("T")-transmission("In");  # get fraction of power absorbed in Silicon layer
SA = getdata("source","area");  # source area, required for final integral

plot(lambda_sim*1e9,Pabs,"wavelength (nm)","normalized power","Absorption");


################################################
# get solar spectrum data


# get solar spectrum as a function of wavelength
lambda= solar(0);  # units of m
#lambda=lambda(length(lambda):-1:1);
ssp    = solar(1);  # units of W/m^2/nm
#ssp=ssp(length(ssp):-1:1);

# trim solar spectrum data to range of simulation data
lmin = find(lambda,min(lambda_sim));
lmax = find(lambda,max(lambda_sim));
lambda = lambda(lmin:lmax);
ssp = ssp(lmin:lmax);

plot(lambda*1e9,ssp,"wavelength (nm)","solar spectrum (W/m^2/m)");


################################################
# Do some analysis using the solar spectrum data

# interpolate simulation data to use the same frequency vector as solar spectrum
Pabs = interp(Pabs,lambda_sim,lambda); 

# Integrate over wavelength to get total power, and 
# multiply by the source area to get the absorbed power in Watts
Pabs_total  = integrate(Pabs  * ssp,1,lambda) * SA;

?"Power absorbed in the simulation volume: " + num2str(Pabs_total) + " Watts";