%=======================================================================
% Radiative parameters (Photosphera) C Aleksander S. Zaitsev | NSU 2000
%=======================================================================

cd E:\USERDAT\Alex\Databs-6s\Sun2000\AtmModel %data files storage

dataset = load('table1a');

%-------------------------------------
n_electr  = 2;
n_ion     = 2;
n_atom    = 4;
k			 = 11606;
I_H       = 13.6056981; 		  %eV
el_charge = 4.8032068*10^10; 	  %SGS
plank		 = 6.6260755*10^(-27); %erg*s
melectr   = 9.1093897*10^(-28); %SGS
stephan   = 5.67051*10^(-5);    %erg cm-2 c-1 K-4
r_sun     = 6.9599*10^10;       %cm
S_sun     = 4*pi*r_sun^2;       %cm2
L_sun     = 3.826*10^33;        %erg s-1
%-------------------------------------

dheight = dataset(:,1);
dtemp   = dataset(:,2);
delectr = dataset(:,3);
dhidr   = dataset(:,4);

alpha = 1./(dhidr./delectr + 1);

lambda  = plank./sqrt(2*pi*melectr*dtemp/k*1.60217733*10^(-12));
tempval = n_electr*n_ion/n_atom*exp(-I_H*k./dtemp)./(delectr + dhidr)./lambda.^3;

saha = tempval.*(sqrt(1+4./tempval)-1)/2;

figure(1);
hold on
  plot(dheight,saha );
  plot(dheight,alpha);
hold off  

delta  = (max(dheight)-min(dheight))/max(size(dheight))*10^5;
p_temp = S_sun*stephan*dtemp.^4;
p_torm = S_sun*1.6*10^(-25)*delectr.^2.*sqrt(dtemp/k)*delta;
p_rec  = S_sun*5  *10^(-24)*delectr.^2./sqrt(dtemp/k)*delta;

figure(2);
hold on
  plot(dheight(1:33),p_temp(1:33));
  plot(dheight(1:33),p_torm(1:33));
  plot(dheight(1:33),p_rec (1:33));
  plot(dheight(1:33),L_sun);
hold off

l_rad  = 2.5*10^37*(dtemp/k).^(7/2)./delectr.^2;
h_edge = (400 - dheight)*10^5;

figure(3);
hold on
  plot(dheight(22:45),l_rad (22:45));
  plot(dheight(22:45),h_edge(22:45));
hold off  