Документ взят из кэша поисковой машины. Адрес оригинального документа : http://crydee.sai.msu.ru/ftproot/users/vab/oscil/old/Eigen_lfix.m Дата изменения: Fri Apr 29 15:13:22 2005 Дата индексирования: Mon Dec 24 05:25:48 2007 Кодировка: Windows-1251

%Загрузка модели 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
%Вычисление других переменных
%
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);
%Вычисление determinant дифференциальных уравнений
lorder = 1;
r_in=0.01; % inner boundary point
r_cut=0.98; % cut point
r_out=ro(1);
%omega in mHz
%list - probe values of frequencies
l_list=length(list);
list_om=[];
for ii=1:l_list,
omega = list(ii);
omega = 2*pi*omega/1e3;
omegao2 = (omega^2)*(Rsun^3)/(Grav*Msun);
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,:));
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*(Grav*Msun)/Rsun^3)/(2*pi)*1e3;
disp([knots r_om r_om-list(ii)]);
list_om=[list_om; r_om];
hold off
end
%save list1g list_om