% This script is for saving of eigen-solution for frequency from list
% 
global go ro 
global lorder lw
global pp_c2 pp_N2 pp_g pp_ro
%
global sol_out sol_in kk_out
global r_in r_cut r_out
mods=load('model-s.txt'); 
%Константы (Книга такая-то)
Grav = 6.670e-8;
Msun = 1.98e33;
Rsun = 6.9599e10;
GMsun = Grav*Msun;
%Передача переменным значений из модели 
r = mods(:,1);
p = mods(:,4);
rho = mods(:,5);
g1 = mods(:,10);
A = mods(:,15);
lq = mods(:,2);
clear mods
% Other variables from model
c2 = g1.*p./rho;
N2 = (A.*exp(lq)*GMsun)./(r.^3);
g = (1./r.^2).*(GMsun*exp(lq));
%Обезразмеривание величин
R3GM = (Rsun^3)/(GMsun);
ro = r/Rsun;
c2o = c2*(Rsun/(GMsun));
N2o = N2*R3GM;
rhoo = rho*(4*pi*Rsun^3/Msun);
go = g*(Rsun*Rsun)/(GMsun);
% interpolational structures
pp_c2 = interp1(ro,c2o,'cubic','pp');
pp_N2 = interp1(ro,N2o,'cubic','pp');
pp_g = interp1(ro,go,'cubic','pp');
pp_ro = interp1(ro,rhoo,'cubic','pp');
%Степень моды
lorder = 1; lw = lorder*(lorder+1);
r_in=0.01; % inner boundary point
r_cut=0.98; % cut point
r_out=ro(1);
%omega in mHz from list - values of frequencies
l_list=length(list);
dd=[];
norm2_erg=zeros(l_list,1);
knots=zeros(l_list,1);
ind_r=find(ro>=r_in);
Xfin = ro(ind_r);Xfin=Xfin(end:-1:1);
Y2_save=zeros(length(Xfin),l_list);
%formation of gaussian profile
gau_pr=exp(-((Xfin-0.8)/0.05).^2);
for ii=1:l_list,
omega = list(ii);
omegao2 = ((2*pi*omega/1e3)^2)*R3GM;
%
dd=deter(omegao2);
[Y_eig1, Y_eig2]=solut_oscil(Xfin);
%
ny1=Y_eig1'.*sqrt(ppval (Xfin,pp_ro))/omegao2;
ny2=Y_eig2'.*sqrt(ppval (Xfin,pp_ro)).*Xfin;
nerg=ny1.^2*lw+ny2.^2;
pp_erg= interp1(Xfin,nerg,'cubic','pp');
norm2_erg(ii)=quad(@ppval,r_in,r_out,[],[],pp_erg);
Y2_save(:,ii)=Y_eig2*sqrt(norm2_erg(ii));
% decomposition guass profile
cnv_gau=sqrt(norm2_erg(ii))*Y_eig2'.*gau_pr.*ppval (Xfin,pp_ro).*(Xfin.^2);
pp_cnv= interp1(Xfin,cnv_gau,'cubic','pp');
am_cnv(ii)=quad(@ppval,r_in,r_out,[],[],pp_cnv);
disp (num2str(norm2_erg(ii)));
%plot(Xfin,ny1,'r',Xfin,ny2,'b',Xfin,nerg,'g');
%
% number of knots or order of mode
knots(ii)=length(find(Y_eig2(1:end-1).*Y_eig2(2:end)<0));
%
disp (['om=' num2str(omega) ' dd=' num2str(dd) ' knots=' int2str(knots(ii))]);
%
end
%plot (knots,norm2_erg);
plot (knots,am_cnv);