% This script is for saving of eigen-solution for given list of freq
%Загрузка модели S
global c2o N2o go ro lorder
%global Xfin_i Yfin_i Xfin_o Yfin_o  % results of current solutions
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;
%Передача переменным значений из модели 
r = mods(:,1);
p = mods(:,4);
rho = mods(:,5);
g1 = mods(:,10);
A = mods(:,15);
lq = mods(:,2);
clear mods
%Степень моды
lorder = 1;
r_in=0.01; % inner boundary point
r_cut=0.98; % cut point
r_out=ro(1);
% Other variables from model
c2 = g1.*p./rho;
N2 = (A.*exp(lq)*Grav*Msun)./(r.^3);
g = (1./r.^2).*(Grav.*(exp(lq)*Msun));
%Обезразмеривание величин
ro = r/Rsun;
c2o = c2*(Rsun/(Grav*Msun));
N2o = N2*((Rsun^3)/(Grav*Msun));
rhoo = rho*(4*pi*(Rsun^3)/Msun);
go = g*(Rsun*Rsun)/(Grav*Msun);
%Степень моды
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=[];Yeig1=zeros(l_list,999);Yeig2=zeros(l_list,999);
for ii=1:l_list,
omega = list(ii);
omega = 2*pi*omega/1e3;
omegao2 = (omega^2)*(Rsun^3)/(Grav*Msun);
%
dd=deter(omegao2);
%
Xfin_i = linspace(r_in,r_cut,950);
Yfin_i = deval(sol_in,Xfin_i);
%
Xfin_o = linspace(r_cut,r_out,50);
Yfin_o = deval(sol_out,Xfin_o);
Yfin_o = Yfin_o*kk_out;
%
Yfin=[Yfin_i Yfin_o(:,2:end)];
Xfin=[Xfin_i Xfin_o(:,2:end)];
plot(Xfin,Yfin);
% number of knots or order of mode
knots=length(find(Yfin_i(2,1:end-1).*Yfin_i(2,2:end)<0));
knots=knots+length(find(Yfin_o(2,1:end-1).*Yfin_o(2,2:end)<0));
disp ([omega dd knots]);
Y_eig1(ii,:)=Yfin(1,:);
Y_eig2(ii,:)=Yfin(2,:);
hold off
end

