function glaze_programm
% set initial variables
clear 

format long

u0=200; %            eruption velocity [m/s]
r0=100; %100;%           vent size [m]
n0=0.03; %           initial gas mass fraction
thetap0 = 1000; %           eruption temperature [K]
z0= 0; %              height of vent above ground [m]
g= 9.81; %           acceleration due to gravity [m/s^2]
omega = 0;    %          time constant for condensation [s]
al = 0.09; %            entrainment factor
l = 0; %              latent heat of water [J/kg] 2.257d06
r_v = 461; %            gas constant for volcanic gas (water vapor) [J/kg/K]
r_d = 287; %            gas constant for dry air [J/kg/K]
rho_l= 1000; %      density of liquid water [kg/m^3]
rho_s = 1000; %           density of pumice [kg/m^3]
c_d = 998; %            specific heat for dry air at constant P [J/kg/K]
c_v = 2000; %           specific heat for volcanic gas at constant P (water vapor) [J/kg/K]
c_s = 920; %            specific heat for for solids at constant P [J/kg/K]
c_l = 4200; %           specific heat for liquid water at constant P [J/kg/K]
p0 = 1e+05; %             pressure at vent level
hstep =100; %           intergration constant for Runge Kutta [m]
%50.d00             height interval for output in file
%2                  1 read from file, 2 US Standart Atmosphere
%0                 with gas thrus =1 without = 0

  % [phi_av,phi_ad, r_ab, c_ab, rho_ad, rho_av, rho_ab] = deal(0);
  
  z_end=8000;
  nstep=ceil(z_end-z0)/hstep;
  
%set some addtitional values
     % heighi =  heighw ;
      ep = r_d/r_v ;
      cdcof = 0.75 ;
      u = u0   ;
      r = r0   ;
      z = z0  ;
      thetap = thetap0  ;
    %  nbh = 0  ;

     [thetaa,p]=usstand_nested(z);

%set initial conditions
      rho_d = p0 /(r_d*thetap);
      phi_d = 0;
      phi_l = 0;

% rho_v = p*n*rho_s / ( r_v*thetap*n*rho_s-p*(1-n) );
      rho_v = p0/(r_v*thetap);  % 
      phi_s = (1 - n0)*rho_v / ( n0*rho_s + (1-n0)*rho_v );
      phi_v = 1 - phi_s - phi_d - phi_l;

% set initial values for Runge Kutta   
    m_d=rho_d*u*r^2*phi_d;%m_d=zeros(nstep,1); 
    m_v= rho_v*u*r^2*phi_v;%m_v=zeros(nstep,1); 
    m_l=rho_l*u*r^2*phi_l;%m_l=zeros(nstep,1); 
    m_s = rho_s*u*r^2*phi_s;
    m = m_d+m_v+m_l+m_s; 
     
  % set initial condition for diff. equations, for indexes of equations see
    %subroutine glaze
    rho_b = rho_d*phi_d + rho_v*phi_v + rho_l*phi_l + rho_s*phi_s;
    c_b = ( m_d*c_d+m_v*c_v+m_l*c_l+m_s*c_s ) / m ;
    erupra = r^2*u*3.14159*rho_b ;
    thermflux = erupra*(thetap-thetaa)*c_b;

    v=rho_b*u^2*r^2;%v=zeros(nstep,1); 
    t= rho_b*u*r^2*c_b*thetap ;%t=zeros(nstep,1); 
    
  for istep=1:nstep    

      y97(1) = m_d;  %m_d(1)=rho_d*u*r^2*phi_d;
      y97(2) = m_v;%m_v(1)= rho_v*u*r^2*phi_v;
      y97(3) = m_l; %m_l(1)=rho_l*u*r^2*phi_l;
      y97(4) = v; %v(1)=rho_b*u^2*r^2;
      y97(5) = t; %t(1)= rho_b*u*r^2*c_b*thetap ;
      %y98(1) = m_d
      %y98(2) = v(istep)
      %y98(3) = t
      
options = odeset('Refine', 1,'RelTol',1e-2,'AbsTol',1e-2); %'MaxStep',0.1);
[zh,y97n] = ode45(@glaze_97_nested,[z z+hstep],y97, options);
   
        m_d = y97n(end,1); %hier setze ich den letzten Eintrag der Berechnung wieder als Startwert für den nächsten Runge-Loop ein
        m_v = y97n(end,2);
        
         if m_v <= 0
             m_v = 0;
         end
%          
        m_l = y97n(end,3);
        
        v = y97n(end,4);
        
        if v <= 0
            v=0;
        end
        
        t = y97n(end,5);
        
        if t <= 0
            t=0;
        end
        
        

% calculate u, r, etc at new height
      m = m_d+ m_v + m_l + m_s;
      if m<=0
          m=0;
      end
      u = v/ m;
      if u <= 0
          u=0;
      end
      c_b = ( m_d*c_d + m_v*c_v + m_l*c_l + m_s*c_s ) / m;
      %rho_b = rho_d*phi_d + rho_v*phi_v + rho_l*phi_l + rho_s*phi_s;
      rho_b = p /(r_v.*t).* (m.^2 .* c_b ./( m_v + ep.*m_d ));
      rn = ( m.^2 ./ ( rho_b.*v) ).^0.5;
      if rn <= 0
          rn=0;
      end
      rcut = (rn-r)./r;
      if rcut <=0
          rcut=0;
      end
      r = rn;
      if r<=0
          r = 0;
      end
      thetap = t ./ ( rho_b.*u.*r.^2.*c_b ); % da stand un anstatt u??????????????????????
      if thetap <=0
          thetap=0;
      end
      phidpphiv = 1 - p./(r_v.*t).*(m.*c_b)./(m_v+ep.*m_d).*( m_l./rho_l+m_s./rho_s );
      phi_l= 1/rho_l.*(p/r_v.*t).*(m_l.*m.*c_b)./(m_v+ep.*m_d);
      phi_s = 1-(phidpphiv)-phi_l ;

      z = z+hstep;
      
%plot erster versuch mit temperatur   

aa=-z:z; %imagesc(aa,r,t)
plot(aa,r)
%xlabel('z')
%ylabel('r')
hold on

% h(1)=subplot(5,1,1);
%    plot(zh,m_d)%  plot(z_plot,y97n(:,1))
%    ylabel('m_d')
%  xlabel('z')
%  hold on
%   h(2)=subplot(5,1,2);
% plot(zh,m_v)%  plot(z_plot,y97n(:,2))
%   ylabel('m_v')
%   xlabel('z')
%   hold on
%  
%   h(3)=subplot(5,1,3);
% plot(zh,m_l)%  plot(z_plot,y97n(:,3))
%   ylabel('m_l')
%   xlabel('z')
%   hold on
%   h(4)=subplot(5,1,4);
% plot(zh,v)%  plot(z_plot,y97n(:,4))
%   ylabel('v')
%   xlabel('z')
%   hold on
%   h(5)=subplot(5,1,5);
% plot(zh,t);% loglog(zh,t); 
%   ylabel('T')
%   xlabel('z')
%   hold on



  end
  
function dy = glaze_97_nested(z,y)
    
    dy=zeros(5,1);

    m1 = y(1)+y(2)+y(3)+m_s; % gesamtmasse in column
    c_b = (y(1)*c_d+y(2)*c_v+y(3)*c_l+m_s*c_s)/m1; % bulk specific heat

    [thetaa,p,wa] = usstand_nested(z);

    phi_av = wa / ( wa + ep );
    phi_ad = 1 - phi_av ;
      
    c_ab = ( c_d + wa*c_v ) / ( 1+wa );
    r_ab = ( r_d + wa*r_v ) / ( 1+wa );
    
    rho_ad = p./( r_v.*thetaa.*(wa+ep) ) ./ phi_ad;
    rho_av = ep*rho_ad;
    rho_ab = rho_av*phi_av + rho_ad*phi_ad;
      
    a1=r_v*y(5)/p;
    b1=y(2)+ep*y(1);
    c1 = b1/(m1^2*c_b);
    d1 = m1/y(4);
    e1=g/y(4);
    f1=a1*b1/c_b;
    g1=b1/(m1*c_b);  %b1/(m1*c_b)
    h1=a1*b1/m1*c_b;
    i1=y(3)/rho_l;
    j1=m_s/rho_s;
% 
    dy(1)=2*al*rho_ad*phi_ad* ((r_v*y(5)/p) *y(4)* ( ( y(2)+ep*y(1)) / (m1^2*c_b) ) )^0.5;
    dy(2)=2*al*rho_av*phi_av*( r_v*y(5)/p *y(4)*(( y(2)+ep*y(1) )/(m1^2*c_b)) )^0.5 - omega*y(2)* m1/y(4);
    dy(3)=omega*y(2)*m1/y(4);
    dy(4)=g/y(4)* (rho_ab* ((r_v*y(5)/p) * ( y(2)+ep*y(1) )/c_b)-m1^2);
    dy(5)= 2*al*rho_ab*c_ab*thetaa*(((r_v*y(5)/p)*(y(4)*(y(2)+ ep*y(1))/(m1^2*c_b)))^0.5)+ l*omega*y(2)*m1/y(4) -rho_ab*g*((r_v*y(5)/p)*((y(2)+ep*y(1))/(m1*c_b))-(m_s/rho_s+y(3)/rho_l)); 
 
    
   % dy(5) = 2*al*rho_ab*c_ab*thetaa * ( (r_v*t/p) * y(4) * (( y(2) + ep*y(1)) /(m1^2*c_b)) )^0.5 + l*omega*y(2)*m1/y(4) -...
   % rho_ab*g*( (r_v*t/p) *  ((y(2) +ep*y(1)) /(m1*c_b)) - (y(3)/rho_l + m_s/rho_s) );
    
    
%     dy(1)=2*al*rho_ad*phi_ad*(a1*y(4)*c1)^0.5;
%     dy(2)=2*al*rho_av*phi_av*(a1*y(4)*c1)^0.5-omega*y(2)*d1;
%     dy(3)=omega*y(2)*d1;
%     dy(4)=e1*(rho_ab*(f1)-m1^2);
%     %dy(5)= 2*al*rho_ab*c_ab*thetaa*(((r_v*y(5)/p)*(y(4)*(y(2)+ ep*y(1))/(m1^2*c_b)))^0.5)+ l*omega*y(2)*m1/y(4) -rho_ab*g*((r_v*y(5)/p)*((y(2)+ep*y(1))/(m1*c_b))-(m_s/rho_s+y(3)/rho_l)); 
%     dy(5)= 2*al*rho_ab*c_ab*thetaa*((a1*(y(4)* c1))^0.5)+ l*omega*y(2)*m1/y(4) -rho_ab*g*(a1*g1-(j1+i1)); 
  
end

% calculate the atmospheric temperature and presssure for a
% US Standart atmosphere

function [te,pr,wa] = usstand_nested(z)

      ta0 = 288.86;
      p0 = 1.01325e+5;
      h1= 11;
      h2= 20;
      xmu1 = 6.5;
      xmu2 = -2.0;

      zkm = z / 1000 ;
      exp1 = g*1000  / (r_d*xmu1);
      exp2 = g*1000 / (r_d*xmu2);

      if  zkm < h1 
        te = ta0-xmu1*zkm;
        pr = p0*((te/ta0)^(exp1));
       elseif ( (zkm >= h1) && (zkm <= h2) )
        te = ta0-xmu1*h1
        pr = p0*((te/ta0)^(exp1))*exp(g*(h1-zkm)*1000.d00/(r_d*te)) % dexp
      elseif  zkm > h2 ; 
        te = ta0-xmu1*h1-xmu2*(zkm-h2)
        c1 = exp(g*(h1-h2)*1000.d00/(r_d*(ta0-xmu1*h1))); %dexp
        c2 = (te/(ta0-xmu1*h1))^(exp2);
        pr = p0*((ta0-xmu1*h1)/ta0)^(exp1)*c1*c2
      end
      
      if zkm > 50
          print 'Plume to high'
      end 
     wa=0;
   

end    
      
end