%
%
%
% Observations from Bokhove 2007/2008 referring to data from Thomas Birner (Birner 2006)
% By Onno Bokhove April 2010
% Use at your own risk; no claims can be made whatsoever; please spy but write your own code.
% In student version parts of this program have been removed (without notice) and the functions Al.m,Fl.m, and bl.m are absent.
% Also, the code can be made to run faster by removing for-loops and using more vector calculations.
% Check everything please.
% 
%
global Ng alp invkap sentrop sB sT Nel rcoeff Mh Mhi dsdzet lam1 lam2 dpsi1 dpsi2 dpsi3 dpsi4 dpsi gg ww psi touw invkap0
%
%
%
cp = 1004.6;
R = 287.04;
sr = R;
theta0 = 290;
Tr = theta0;
theta10 = 320;
theta20 = 500;
s0s = sr+cp*log(theta0/Tr)
s10s = sr+cp*log(theta10/Tr)
s20s = sr+cp*log(theta20/Tr)
%
alp = sr/R
kappa = R/cp
invkap = 1/kappa;
s0 = s0s/sr   % 1
s10 = s10s/sr % about 1.34
s20 = s20s/sr % about 2.91
touw = (theta10-theta0)/(s10-s0)/Tr;
sB = s0;
sT = s10;
pT = (1.0+touw*(sT-1))^(invkap)*exp(-alp*(sT-1.0));
Ng = 5;
gg = [-sqrt(5+2*sqrt(10/7))/3 -sqrt(5-2*sqrt(10/7))/3 0 sqrt(5-2*sqrt(10/7))/3 sqrt(5+2*sqrt(10/7))/3]; % Gaussian integration points in reference coordinate interval -1<zeta<1 (from web)
ww = [(322-13*sqrt(70))/900 (322+13*sqrt(70))/900  128/225 (322+13*sqrt(70))/900 (322-13*sqrt(70))/900 ]; % Gaussian weights
lam1 = 0.5*(1-gg);
lam2 = 0.5*(1+gg);
dpsi1 = -3.0*lam1.*lam2;
dpsi2 = 3*lam1.*lam2; 
dpsi3 = lam1.*(3*lam1-2.0); 
dpsi4 = lam2.*(3*lam2-2.0); 
dpsi = [dpsi1 dpsi2 dpsi3 dpsi4];
psi = [(lam1.^2).*(2.0*lam2+1) (lam2.^2).*(2.0*lam1+1) 2.0*lam1.^2.*lam2 -2.0*lam2.^2.*lam1];
rcoeff = [sB sT pT Ng gg ww];
%
%
%
Nel = 100; % number of 1D elements
Nno = Nel+1; % number of nodes in 1D
Mh = zeros(2*Nno,1);
Mhi = zeros(2*Nno,1);
sentrop = sB:(sT-sB)/Nel:sT;
dsdzet = sentrop(2:Nno)-sentrop(1:Nno-1);
e = ones(2*Nno,1);
Ag = spdiags(0*[e e e e e e e], -3:3, Nno, Nno); % diagonal matrix to be used
spy(Ag)
clear e;  
invkap0 = invkap;
% 
% Construction of global RHS vector bg
% 
k = 1;
l = 2*k-1;
bg(l) = bl(1,k);
l = 2*k;
bg(l) = bl(3,k);
for k = 2:Nel
  l = 2*k-1;
  bg(l) = bl(2,k-1)+bl(1,k);
  l = 2*k;
  bg(l) = bl(4,k-1)+bl(3,k);
end
k = Nel+1;
l = 2*k-1;
bg(l) = bl(2,k-1)+pT;
l = 2*k;
bg(l) = bl(4,k-1);
%
%
% Construction of global matrix
%
%
% Newton/Picard iteration starts
%	(Picard did not seem to work; start up with linear case: invkap0=1 then rest; first test with invkap=1 only)
%
%
nfix = 1; invkap0 = 1;
figure(12); clf;
nit = 1; % iterate counter
Mh = Mhi+10^6;
while (nit < 2 | max(abs(Mhi-Mh))>1.e-9 ) % quit test
 % 
 Mh = Mhi;
 nit = nit+1
 switch nfix
 case 2
  k = 1;
  l = 2*k-1;
  F0(l) = Fl(1,k)+Mhi(1)^invkap-bg(l); % F0=global vector estimate of weak form using prior solution
  l = 2*k;
  % REMOVED
  %
 end
 %
 k = 1;
 l = 2*k-1;
  Ag(l,2*k-1) = Al(1,1,k,nit)+Mhi(1)^(invkap0-1); % Ag=global matrix used in Picard iteraion=nearly matrix used in Newton itearion
  Ag(l,2*k)   = Al(1,3,k,nit); % Al = local element matrix for Ag
  Ag(l,2*k+1) = Al(1,2,k,nit);
  Ag(l,2*k+2) = Al(1,4,k,nit);
 l = 2*k;
 % REMOVED
 %
 %
 %
 switch nfix
 case 1 % fixed point iteration (does not seem to work/converge)
   Mhi = Ag\bg';
   nfix = 2;
   invkap0 = invkap;
 case 2 % Newton iteration
   Mhi = Mhi-kappa*(Ag\F0'); % Check!
 end
 for k = 1:Nel+1
   Moh(k)   = Mhi(2*k-1); % M-values
   dMhds(k) = Mhi(2*k);   % derivatives; note confusion in reader: dM/ds should be one at some nodes, not dM/dzeta!
 end
 %
 figure(12); % Plotting of iterates
 %
 subplot(2,1,1);
 % plot(sentrop,Moh,'-r'); hold on;
 plot(sentrop,Moh,'.b'); hold on;
 plot(sentrop,1+(1-touw)*(sentrop-1)+0.5*touw*(sentrop.^2-1),'-k'); hold on;
 subplot(2,1,2);
 % plot(sentrop,dMhds,'-b'); hold on;
 plot(sentrop,dMhds,'or'); hold on;
 plot(sentrop,(1-touw+touw*sentrop),'-k'); hold on;
 pause(0.5);
 invkap0 = invkap;
end