function tide4func_fin
%'final' version, can use any .GYRE data
%use global variable, avoid readin data in every function calls

%fix the eigenfreq to the forcing frequency, then get eigenfunctions
%following GYRE paper

%--------------------------readin the variables---------------------
L=2.0;

Rsun=6.955D10;
Msun=1.9891D33;
Grav=6.67428D-8;%the value used in GYRE in file
pai=3.1415926535898D0;
P43=4.0D0*pai/3.0D0;

arr1=dlmread('profile181.data.GYRE');
%already deleted the first line

%the general parameter:Mstar,Rstar in the deleted first line in profile181.data.GYRE
%Mstar=5.0*Msun;
Mstar=0.994600000000E34; %the same,also 5Msun
Rstar=0.182074753876E12;
M2=5.0*Msun;

nz=895;
%k r w Lr p T rho Delta N2 gamma1 Delta_ad del kappa kappa_T kappa_rho eps eps_T eps_rho Omega_rot
%1 2 3 4  5 6  7    8    9  10     11       12  13     14      15       16   17    18       19	
k=arr1(1:nz,1);
r=arr1(1:nz,2);
w=arr1(1:nz,3);%Mr/(Mstar-Mr)
Lr=arr1(1:nz,4);
p=arr1(1:nz,5);
T=arr1(1:nz,6);
rho=arr1(1:nz,7);
Delta=arr1(1:nz,8);
N2=arr1(1:nz,9);
gamma1=arr1(1:nz,10);
Delta_ad=arr1(1:nz,11);
del=arr1(1:nz,12);
%kappa=arr1(1:nz,13);
%kappa_T=arr1(1:nz,14);
%kappa_rho=arr1(1:nz,15);
%eps=arr1(1:nz,16);
%eps_T=arr1(1:nz,17);
%eps_rho=arr1(1:nz,18);
%omega_rot=arr1(1:nz,19);

%note the ./ .* .^ operator needed in matlab
x=r/Rstar;
Mr=w./(1.0+w)*Mstar;
g=Grav.*Mr./r.^2;

c2=gamma1.*p./rho;
S2_L=L*(L+1.0).*c2./r.^2;
Astar=r./g.*N2;
U=4.0*pai*rho.*r.^3./Mr;
Vg=g.*r./c2;
V=Vg.*gamma1;
%c1=x^3*(1.0+w)/w
%omg2=sig^2*Rstar^3/(Grav*Mstar); %dimless freq^2
%-------------------------
y1=g./r;
y2=L*(L+1.0).*g./r./S2_L;
%y3=Astar
%y4=U
%y5=1.0./(1.0+V)
%----------
%V=1.0./y5-1.0; %note,need element divide './'
%U=y4;
%gamma1=V./y2;
%Vg=V./gamma1;
%Astar=y3;

rho_less=rho/(Mstar/Rstar^3);
g_less=g/(Grav*Mstar/Rstar^2);


clmk=-0.5*0.7335169487139; %eq 1.37 in Willems_PhD
%for l=2,m=0,k=1 from hansen4a.pro
%clmk=c201=1*Plm(0)*X_lmk=-0.5*0.7335169487139
%X_lmk in Burkart is the same as the X_-(L+1)_-m in Willems_PhD
%Plm(0)=1/sqrt(pi)*Gamma(1.5)/Gamma(2)*cos(pi)
%-------------------
realx=log(r/Rstar);
dlmwrite('realx.txt',realx,' ');

global Rsun,Msun,Grav,pai,P43,Mstar,Rstar,M2,nz,k,r,w,Lr,p,T,rho,Delta,N2,gamma1,Delta_ad,x,Mr,g,c2,S2_L,Astar,U,Vg,V,y1,y2,rho_less,g_less,clmk,realx 

%--------------------------run the code ------------------------------change this part to run the code
%this pro uses a 5Msun MS star R=2.62Rsun,
%calculate the eigenfunction using(Dziem outerBC),around the g9 mode, Period=44448.0, 
%see mesa/5.0M/gmode/DziemBC/tide5D.pro
% P=48726.537sec(g10 mode)  ,P=44448.888sec(g9 mode) 
% the forcing period vs eps_r(1) diagram, see plotepsr.pro; similar to Fig5 in CAFein paper
p0=44448.0;
pend=44450.0;
%p0=48500.0;
%pend=48900.0;
np=2;
guessPer=linspace(p0,pend,np);
%freq=1/Per
%angular freq=2*pai/Per
%dimless angular freq=wT=2*pai/Per/sqrt(Grav*Msar/Rstar^3)
for i=1:np
temp=guessPer(i);
re=detM( temp );
comb=[temp-44400.0,re];
disp(comb);
%fprintf('per =%25.19f.\n',temp)
%fprintf('per =%25.19f.\n',re)
end

%-------------------------------------------------------------------
function result = detM(guessPeriod)
%return the abs(eps_r) at surface
L=2.0;
x0=realx(1); %center
xend=realx(nz); %surface

%get a from P,and M1,M2
a=(Grav*(Mstar+M2)*guessPeriod^2/(4.0*pai^2))^(1.0/3.0);
epsT=(Rstar/a)^3*(M2/Mstar);


%return the abs(eps_r) at surface
sig=2.0*pai/guessPeriod; %freq
c1omg2=sig^2./y1;
omg2=sig^2/(Grav*Mstar/Rstar^3);%dimless angular freq

%rho_0=interp1(realx,rho,x0,'spline');
%rho_0=rho(1);

%better to move all these quantities outside the function and use global
%global r,g,rho,y1,y2,y3,y4,y5,v,u,gamma1,Vg,Astar,Rsun,L,realx,rho_0,Grav,P43,pai;


%BC at center, CAFein paper eq(15,16,61)
Bin=zeros(3,6);
Bin(1,1)=c1omg2(1);
Bin(1,2)=-1.0*L;
%Bin(1,3)=0.0;
%Bin(1,4)=0.0;
%Bin(1,4)=0.0;
%Bin(1,4)=0.0;

%Bin(2,1)=0.0;
%Bin(2,2)=0.0;
Bin(2,3)=L;
Bin(2,4)=-1.0;
%Bin(2,5)=0.0;
%Bin(2,6)=0.0;

%Bin(3,1)=0.0
%Bin(3,2)=0.0
%Bin(3,3)=0.0
%Bin(3,4)=0.0
Bin(3,5)=1.0;
Bin(3,6)=-1.0;
%print,'Bin= ',Bin
%----------------------
Vsurf=interp1(realx,V,xend,'spline');
%BC at surf, CAFein paper eq(18,59,61)
Bout=zeros(3,6);
Bout(1,1)=1.0+(L*(L+1.0)/omg2-4.0-omg2)/Vsurf;
Bout(1,2)=-1.0;
Bout(1,3)=1.0+(L*(L+1.0)/omg2-L-1.0)/Vsurf;
%Bout(1,4)=0.0;
%Bout(1,5)=0.0;
%Bout(1,6)=0.0;

%CAFein paper eq(59)
Bout(2,1)=4.0*pai*rho_less(nz)/g_less(nz);
Bout(2,2)=0.0;
Bout(2,3)=L+1.0;
Bout(2,4)=1.0;
Bout(2,5)=0.0;
Bout(2,6)=epsT*(2.0*L+1.0)*clmk/g_less(nz);

%Bout(3,1)=0.0
%Bout(3,2)=0.0;
%Bout(3,3)=0.0;
%Bout(3,4)=0.0;
Bout(3,5)=1.0;
Bout(3,6)=-1.0;

%print,'Bout=',Bout
%disp(rho_less(nz));
%disp(g_less(nz));
fprintf('per =%25.19f.\n',rho_less(nz))
fprintf('per =%25.19f.\n',g_less(nz))
fprintf('per =%25.19f.\n',clmk)
%--------------------only 1->N-1 are needed
%coefficient of y123478
T_arr=zeros(6,6,nz-1);
%[ (Vg-3), 						( L*(L+1.0D0)*y1/sig^2-Vg)  	 Vg    		 		 0		0		0
% (sig^2/y1-Astar)  				 (1.-U+Astar)  	    	-Astar      		 0		0		0
% 				0 										0							   (1.0D0-U)    		 1		0		0
%			Astar*U     						 Vg*U    					 L*(L+1.0D0)-Vg*U  -U 	0		0 
%    0													0											0						 0    0		0
%%   0													0											0						 0    0		0]
diag6=[1 1 1 1 1 1];
I66=diag(diag6);

%only need 1-> N-1 for A and T
for i=1:nz-1  
dr=realx(i+1)-realx(i);
xmid=realx(i)+0.5*dr;
Ai=zeros(6,6);

Vg_half=interp1(realx,Vg,xmid,'spline');
y1_half=interp1(realx,y1,xmid,'spline');
Astar_half=interp1(realx,Astar,xmid,'spline');
U_half=interp1(realx,U,xmid,'spline');

Ai(1,1)=Vg_half-3.0;
Ai(1,2)=( L*(L+1.0)*y1_half/sig^2-Vg_half);
Ai(1,3)=Vg_half;
Ai(1,4)=0.0;
Ai(1,5)=0.0;
Ai(1,6)=0.0;

Ai(2,1)=(sig^2/y1_half-Astar_half);
Ai(2,2)=(1.0-U_half+Astar_half);
Ai(2,3)=-Astar_half;
Ai(2,4)=0.0;
Ai(2,5)=0.0;
Ai(2,6)=0.0;

Ai(3,1)=0.;
Ai(3,2)=0.;
Ai(3,3)=1.0-U_half;
Ai(3,4)=1.0;
Ai(3,5)=0.0;
Ai(3,6)=0.0;

Ai(4,1)=Astar_half*U_half;
Ai(4,2)=Vg_half*U_half;
Ai(4,3)=L*(L+1.0)-Vg_half*U_half;
Ai(4,4)=-U_half;
Ai(4,5)=0.0;
Ai(4,6)=0.0;

Ai(5,1)=0.0;
%Ai(5,2)=0.0;
%Ai(5,3)=0.0;
%Ai(5,4)=0.0;
%Ai(5,5)=0.0;
%Ai(5,6)=0.0;

Ai(6,1)=0.0;
%Ai(6,2)=0.0;
%Ai(6,3)=0.0;
%Ai(6,4)=0.0;
%Ai(6,5)=0.0;
%Ai(6,6)=0.0;
%expm(Ai.*dr)
%Ai_arr(:,:,i)=Ai;
%T_arr(:,:,i)=expm(Ai.*dr);
T_arr(:,:,i)=expm_new(Ai.*dr);
end

nm=6*nz;
%nz=895, nm=6*895=5370

MM=zeros(nm,nm);
MM(1:3,1:6)=Bin;
MM(nm-2:nm,nm-5:nm)=Bout;
for i=0:nz-2
 MM(4+i*6:9+i*6,1+i*6:6+i*6)=-T_arr(:,:,i+1);
 MM(4+i*6:9+i*6,7+i*6:12+i*6)=I66;
end

%result=det(M);
con=cond(MM); %this step use svd, so slow
%1.8909e+20

MMb = balance(MM);%,'noperm')
con2=cond(MMb)
%1.6213e+17

%it seems that use the MMb (MM after balance) give similar eigenvectors
[LL,UU] = lu(MMb);
%note, used to use [L,U]=lu(M), so L and U are misconfused with L=2,U=coefficient

%find index of zero elements, should be the right corner one
%for i=1:5370
 %temp2=UU(i,i);
 %if(abs(temp2) < 0.00001)
    %disp(i)
 		%fprintf('Uii =%24.19f.\n',temp2);
 %end		
%end 

jp=nm; %zero elements index,should be the last one,so =nm
Uprim=UU(1:jp-1,1:jp-1);
uuprim=UU(1:jp-1,jp) ;

sol=backsub(Uprim,uuprim);
%[n,m]=size(sol);

%add 1.0 to the end of sol
yall=ones(nm,1);
yall(1:jp-1,1)=sol;

%yall=null(MM); null command to find solution to Ax=0, based on svd
%seems to give results with a sign difference

%eigenfunctions
y123478=zeros(6,nz);
for i=1:nz
	y123478(1,i)=yall(1+6*(i-1));
	y123478(2,i)=yall(2+6*(i-1));
	y123478(3,i)=yall(3+6*(i-1));
	y123478(4,i)=yall(4+6*(i-1));
	y123478(5,i)=yall(5+6*(i-1));
	y123478(6,i)=yall(6+6*(i-1));
end	
y1eig=y123478(1,:);
y2eig=y123478(2,:);
y3eig=y123478(3,:);
y4eig=y123478(4,:);
y7eig=y123478(5,:);
y8eig=y123478(6,:);

epsr_surf=abs(y1eig(nz));
result=epsr_surf;
%disp(epsr_surf)

%observe how eigenfunc change around the free g9 mode, wT=0.42625;  /DziemBC/tide5D.pro
figure
plot(realx,y1eig,'.r')
%hold on

figure
plot(realx,y2eig,'*b')

%figure
%plot(realx,y3eig,'.g')

%figure
%plot(realx,y4eig,'*r')

y123478x=zeros(7,nz);
y123478x(1:6,:)=y123478;
y123478x(7,:)=realx;
num=num2str(guessPeriod);
outname=strcat('y123478x_',num,'.txt');
%dlmwrite(outname,transpose(y123478x),' ');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x=backsub(A,y);
% backsub: solves upper triangular system A*x = y using backwards
%substitution
% usage: x=backsub(A,y);
%
% arguments:
% A (nxn) - upper triangular matrix
% y (nx1) - vector
%
% x (nx1) - solution to A*x = y

% check to see that A is upper triangular
if any(any(tril(A,-1)))
   error('A must be upper triangular')
end

% check dimensions of A and y
[n,m]=size(A);
if n~=m
   error('A must be square');
end
if any([n 1]~=size(y))
   error('y must be a vector with the same # of rows as A');
end

% perform backwards substitution
x=zeros(n,1);
x(n)=y(n)/A(n,n);
for i=n-1:-1:1
   x(i)=(y(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);
end

end
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX



end										