Документ взят из кэша поисковой машины. Адрес оригинального документа : http://crydee.sai.msu.ru/ftproot/users/vab/oscil/Eigen_lfix.m Дата изменения: Mon May 2 16:29:00 2005 Дата индексирования: Tue Oct 2 16:24:56 2012 Кодировка: Windows-1251

% searching eigensolutions from probe values in list
global go ro % used in deter for boundary condition
global lorder lw
global pp_c2 pp_N2 pp_g pp_ro % interpolational structures of models
%
global sol_out sol_in kk_out % result of eigen solution
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); %gamma_1
A = mods(:,15);
lq = mods(:,2); %lg m_r
clear mods
%Вычисление дополнительных переменных
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^2)/(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');
%Вычисление determinant дифференциальных уравнений
lorder = 0; lw = lorder*(lorder+1);
r_in=0.01; % inner boundary point
r_cut=0.98; % cut point
r_out=ro(1);
%
%list - probe values of frequencies
l_list=length(list);
list_om=[];
for ii=1:l_list,
omega = list(ii);
omegao2 = (2*pi*omega/1e3)^2*R3GM;
% searching for zeros of determinant
options=optimset('TolX',1e-5);
om=fzero('deter',omegao2,options);
%
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;
%
plot(Xfin_i, Yfin_i(1,:));
hold on
plot (Xfin_i, Yfin_i(2,:));
plot(Xfin_o, Yfin_o(1,:));
plot (Xfin_o, Yfin_o(2,:));
% a number of zeros in solution
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));
r_om=sqrt(om/R3GM)/(2*pi)*1e3; % return from dimless to mHz
disp([knots r_om r_om-list(ii)]);
list_om=[list_om; r_om];
hold off
end
%save list1g list_om