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

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

dataset = load('table1a'   );
en      = load('e_norm.dat');

%-------------------------------------
n_electr  = 2;
n_ion     = 2;
n_atom    = 4;
k			 = 11606;
I_H       = 13.6056981; 		  %eV
el_charge = 4.8032068*10^10; 	  %SGS
el_si     = 1.60217733*10^-(19);%Kl
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
c         = 2.99792458*10^10;   %cm s-1
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)*10^5/r_sun; %cm
dtemp   = dataset(:,2)/k; %K
n_h     = dataset(:,4)+dataset(:,3); %cm-3

%MESHes:
rNUM = 100;
hNUM = 100;
hMIN = 0  *10^5/r_sun;
hMAX = 300*10^5/r_sun;

r  = 0:(1/rNUM):1;
h  = hMIN:(hMAX-hMIN)/hNUM:hMAX;

for i = 1:1:rNUM
  hl   (i,:) = 1 - sqrt(1+h.^2-2*h.*sqrt(1+r(i).^2));
  temp (i,:) = spline(dheight,dtemp,hl(i,:));
  ntemp(i,:) = spline(dheight,n_h  ,hl(i,:)); 
  alpha(i,:) = ntemp(i,:).*temp(i,:).^3.*exp(-I_H./temp(i,:)).*(1-exp(-I_H./temp(i,:)));
  b    (i,:) = temp(i,:).^4;
end

I = sum(b./alpha,1)';
I = I/max(max(I));

%figure(1)
%mesh(alpha)

%figure(2)
%mesh(hl)

%figure(1)
%hold on
%  plot(r(1:1:max(size(I))),0.970*sqrt(1-I))
%  plot((1 - en(30:600,1))*0.60/1.65 + 0.652,en(30:600,2)/1.05)
%hold off  

xx = 0.65:0.01:1;
r1 = spline(r(1:1:max(size(I))),0.970*sqrt(1-I),xx);
r2 = spline((1 - en(30:600,1))*0.60/1.65 + 0.652,en(30:600,2)/1.05,xx);
hi = sum((1-r1./r2).^2)/max(size(r1))

%-------------------------------------
cd E:\USERDAT\Alex\Databs-6s\Sun2000\