function glaze%(u0,r0,n0,thetap0)
% set initial variables
%clear 

format long

 u0=300;        %eruption velocity [m/s]
 r0=100; %100;%           vent size [m]
 n0=0.02; %           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 =500; %           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

  
  z_end=40000;
  nstep=ceil(z_end-z0)/hstep;
  r_vec=zeros(nstep,1);
%set some addtitional values

     % heighi =  heighw ;
      ep = r_d/r_v ;%0.622
      cdcof = 0.75 ;
      u = u0   ;
      r = r0   ;
      z = z0  ;
      thetap = thetap0  ;
 

     [thetaa,p]=usstand_nested(z);

%set initial conditions
      rho_d = p0 /(r_d*thetap);
      phi_d = 0;
      phi_l = 0;
      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_v= rho_v*u*r^2*phi_v;
    m_l=rho_l*u*r^2*phi_l;
    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;
    rho_b0=rho_b;
    m_v0=m_v;
    c_b = ( m_d*c_d+m_v*c_v+m_l*c_l+m_s*c_s ) / m ;
    
   if m <=0
       c_b=0;
   end
    
    erupra = r^2*u*3.14159*rho_b ;
    thermflux = erupra*(thetap-thetaa)*c_b;

    v=rho_b*u^2*r^2; 
    t= rho_b*u*r^2*c_b*thetap ;
    
igasthrust = 1   ; 

  
      y97(1) = m_d;    
      y97(2) = m_v;    
      y97(3) = m_l;    
      y97(4) = v;              
      y97(5) = t;      

      y98(1) = m_d;
      y98(2) = v;       
      y98(3) = t;      
                    
       
                
  for istep=1:nstep    

      
options = odeset('Refine', 1,'RelTol',1e-1,'AbsTol',1e-1); %'MaxStep',0.1);
[zh,y97n] = ode45(@glaze_97_nested,[z z+hstep],y97,options);
  
options = odeset('Refine', 1,'RelTol',1e-1,'AbsTol',1e-1); %'MaxStep',0.1);
[zh,y98n] = ode45(@glaze_98_nested,[z z+hstep],y98, options);


    if igasthrust ==1
        
               
        m_d  = y98n(end,1); %hier setze ich den letzten Eintrag der Berechnung wieder als Startwert für den nächsten Runge-Loop ein
        v    = y98n(end,2);        
        t    = y98n(end,3);         
                    
   if y97n(4)- y98n(2)>0 && y97n(5)- y98n(3)>0  
        igasthrust =0;
   end          
                    
    else

        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);              
        m_l = y97n(end,3); 
        v   = y97n(end,4);         
        t   = y97n(end,5);                
    end
    
  % calculate u, r, etc at new height

         m = m_d+ m_v + m_l + m_s;   
         u = v/ m;
                    if u <= 5
                    u=0;
                    end  
                    
                                                                
        c_b = ( m_d*c_d + m_v*c_v + m_l*c_l + m_s*c_s ) / m;   
        
        if m<=0;
            c_b=0;
        end
        
        rho_b = p /(r_v*t)* (m^2 * c_b /( m_v + ep*m_d ));
        if t<=0
            rho_b=0;
        end
        rn = ( m^2 / ( rho_b*v) )^0.5;
        if rho_b<=0
            rn=0;
        end
        if v<=0;
            rn=0;
        end
        rcut = (rn-r)/r;
        if r<=0;
            rcut=0;
        end
        r = rn;
        r_vec(istep) = r;
        thetap = t / ( rho_b*u*r^2*c_b ); 
        if rho_b<=0;
            thetap=0;
        end
        if u<=5;
            thetap=0;
        end
        if r<=0;
            thetap=0;
        end
        if c_b<=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 );   
        if r_v <=0;
            phidpphiv =0 ;
        end
        if t <=0;
            phidpphiv =0 ;
        end
     
            
        
        phi_l= 1/rho_l*(p/r_v*t)*(m_l*m*c_b)/(m_v+ep*m_d);
        if r_v<=0;
            phi_l=0;
        end
       if t<=0;
            phi_l=0;
        end 
        phi_s = 1-(phidpphiv)-phi_l ;

        if igasthrust ==1
    
    y98(1)=y98n(end,1);
    y98(2)=y98n(end,2);
    y98(3)=y98n(end,3);
    y97(1)=y98n(end,1);
    y97(2)= m_v;
    y97(3)= 0;
    y97(4)=y98n(end,2);
    y97(5)=y98n(end,3); 
    
else
    
       y97(1)=y97n(end,1);
       y97(2)=y97n(end,2);
       y97(3)=y97n(end,3);
       y97(4)=y97n(end,4);
       y97(5)=y97n(end,5); 
       
        end

 
       z = z+hstep;
            
  end

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

 
hold on

  z_vec = (hstep:hstep:z)';
%   plot (r_vec/10e4,z_vec/10e4,'r-')
%   plot (-r_vec*1/(10e4),z_vec*1/(10e4),'r-')
%   axis ([-3 3 -.5 4]);
% 
  plot (r_vec,z_vec,'r-')
  plot (-r_vec,z_vec,'r-')
  xlabel('Breite (km)')
  ylabel('Höhe über Kraterradius (km)')
  
  hold all
  axis ([-30000 30000 -5000 40000]);

%  index=find(isfinite(r_vec));
  %t_vec=zeros(nstep,1);
 % t_vec(istep)=t;
  % z_temp=z_vec(index);
  % r_temp=r_vec(index);
  %a=length(r_temp);
  %t_temp=t_vec(1:a);
 % imagesc(z_temp,r_temp,t_temp);
        
  z2 = [z_vec; z_vec(end:-1:1)];
  r2= [r_vec;-r_vec(end:-1:1)];
%   z2 = [z_temp; z_temp(end:-1:0.0001)];
%   r2 = [r_temp; -r_temp(end:-1:0.0001)];
  patch(r2,z2,ones(size(r2)),'FaceColor','r');

 
h_vulkan = 5000;

hold on

x =[-h_vulkan              h_vulkan       -h_vulkan           h_vulkan ;          
    -h_vulkan              h_vulkan       -200                200;
     h_vulkan             -h_vulkan       -200                200];

y =[-h_vulkan            -h_vulkan        -h_vulkan+ 2200       -h_vulkan+ 2200;
    -h_vulkan+ 2200      -h_vulkan+ 2200     0                     0;
    -h_vulkan            -h_vulkan+ 2200  -h_vulkan+ 2200       -h_vulkan+ 2200];


% 
 z = [1 1 1 1;
      1 1 1 1;
      1 1 1 1];
  
% 
% z = [0.0001 0.0001 0.0001 0.0001 ;
%      0.0001 0.0001 0.0001 0.0001 ;
%      0.0001 0.0001 0.0001 0.0001 ];
    
patch(x,y,z,'FaceColor','black')



hold on
h_vulkan = 5000;
  

x =[-h_vulkan               h_vulkan        -h_vulkan           h_vulkan ;          
    -h_vulkan               h_vulkan        -200                200;
     h_vulkan              -h_vulkan        -200                200];

y =[-h_vulkan             -h_vulkan       -h_vulkan+2200      -h_vulkan+2200;
    -h_vulkan+2200        -h_vulkan+2200       0                       0;
    -h_vulkan             -h_vulkan+2200  -h_vulkan+2200      -h_vulkan+2200];


 z = [1 1 1 1;
      1 1 1 1;
      1 1 1 1];
  
% z = [0.0001 0.0001 0.0001 0.0001 ;
%      0.0001 0.0001 0.0001 0.0001  ;
%      0.0001 0.0001 0.0001 0.0001 ];
    

patch(x,y,z,'FaceColor','black')

hold on
x = [-200                  -200;
     +200                   200;
     -200                   200];
  
y = [ -h_vulkan+2200         0;
      -h_vulkan+2200        -h_vulkan+2200;
         0                    0];
 
%  z= [0.0001 0.0001;
%      0.0001 0.0001;
%      0.0001 0.0001];

  z = [1 1 ;
      1 1 ;
      1 1 ];
 patch(x,y,z,'FaceColor','red')
 


 
 
if v<=0
    v=0;
end
 
if isnan(v)||isinf(v) 
    disp('break')
     return
end


 if isnan(u)||isinf(u) 
    disp('break')
     return
 end
 

   if u == 0
   disp('break')
    return
   end
   
  if r<=0
    r=0;
  end

 if isnan(r)|| isinf(r)
     disp('break')
     return
 end
 
 
 
  if isnan(rn) || isinf(rn)
     disp('break')    
     return
  end
  
   if m<=0
    m=0;
   end

if isnan(m) || isinf(m)
     disp('break')
     return
end
 
  
  % -----------------  ohne Gasschub  ------------------------------
  
  
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
    if m1<=0;
        c_b=0;
    end

    [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))*(1/(wa+ep) ) )/ phi_ad;
    if r_v<=0;
        rho_ad=0;
    end
    if phi_ad<=0;
        rho_ad=0;
    end
    rho_av = ep*rho_ad; %( (p/( r_v*thetaa))*(wa/(wa+ep)) )/ phi_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);
    m11=m1^2;
    c1 = b1/(m11*c_b);    
    if c_b<=0;
        c1=0;
    end
    d1 = m1/y(4);   
    if y(4)<=0;
        d1=0;
    end
    e1=g/y(4); 
    if y(4)<=0;
        e1=0;
    end
    f1=a1*b1/c_b;
    if c_b<=0;
        f1=0;
    end
    g1=b1/(m1*c_b);
    if m1<=0;
        g1=0;
    end
    if c_b<=0;
        g1=0;
    end
    h11=(a1*b1)/c_b; 
    if c_b<=0;
        h11=0;
    end
    i1=y(3)/rho_l;
    if rho_l<=0;
        i1=0;
    end
    j1=m_s/rho_s;
    if rho_s<=0;
        j1=0;
    end
   

    dy(1)=2*al*rho_ad*phi_ad*(sqrt (a1 *y(4)* c1));
    dy(2)=2*al*rho_av*phi_av*(sqrt(a1 *y(4)*c1)) - omega*y(2)* d1;
    dy(3)=omega*y(2)*d1;   
    dy(4)=e1*(rho_ab* h11-m11);
    dy(5)= 2*al*rho_ab*c_ab*thetaa*(sqrt(a1*y(4)*c1))+ l*omega*y(2)*d1 -rho_ab*g*((a1*g1)-(i1+j1)); 
       
    
   
   
end 
 
   %  --------------   mit gasschubregion  ----------------------
   
   function dy = glaze_98_nested(z,y)
       
    dy=zeros(3,1);

    m2 = y(1)+m_v+m_l+m_s;
   
  
    c_b = (y(1)*c_d + m_v*c_v + m_l*c_l + m_s*c_s)/m2; 

    [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))*(1/(wa+ep) ) )/ phi_ad;
    if r_v<=0;
       rho_ad=0;
    end
    if thetaa<=0;
        rho_ad=0;
    end
    rho_av = ep*rho_ad;%( (p/( r_v*thetaa))*(wa/(wa+ep)) )/ phi_av; % = ep*rho_ad;
    rho_ab = rho_av*phi_av + rho_ad*phi_ad;

    dy(1)= 0.125*(rho_ab*y(2))^0.5;  
    
    dy(2)= g/y(2)*( (rho_ab*r_v*y(3)/p*...
           (m_v+ep*y(1))/((y(1)*c_d+m_v*c_v+m_l*c_l+m_s*c_s)/m2)-m2^2) );
       
    dy(3) =0.125*( rho_ab*y(2) )^0.5 * c_ab * thetaa - rho_ab*g*...
        ( ( (r_v*y(3) )/p) *(m_v+ep*y(1)) /(m2*c_b)-m_s/rho_s );
     
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;
        t1=te/ta0;
        pr = p0*((t1)^(exp1));
       elseif ( (zkm >= h1) && (zkm <= h2) )
        te = ta0-xmu1*h1;
        t2=te/ta0;
        ab1=1000/(r_d*te);
        pr = p0*(t2^(exp1))*exp(g*(h1-zkm)*ab1); % dexp
      elseif  zkm > h2 ; 
        te = ta0-xmu1*h1-xmu2*(zkm-h2);
        ab2=1000/(r_d*(ta0-xmu1*h1));
        c1 = exp(g*(h1-h2)*ab2); %dexp
        ab3=te/(ta0-xmu1*h1);
        c2 = (ab3)^(exp2);
        ab_1=(ta0-xmu1*h1)/ta0;
        pr = p0*(ab_1)^(exp1)*c1*c2;
      end
      
      if zkm > 50
          print 'Plume to high'
      end 
     wa=0;
   

end    
      
end