function glaze(u0,r0,n0,thetap0)
% set initial variables
%clear 
 
format long
% % % 
%   u0=350;        %eruption velocity [m/s]
%   r0=150; %100;%           vent size [m]
%   n0=0.045; %           initial gas mass fraction
%   thetap0 = 1200; %           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);
  t_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 ;
 
        t_vec(istep)=t;
        
        
      
        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
 


%[0.885515873015873 0.127801621363853 0.0158730158730159 0.400095374344301]);
% % %cb=colorbar
 hold all

backgroundImage = importdata('vulkan_Bild.jpg');
%select the axes
%axes(handles.Image);
%place image onto the axes
image([-40 40 ],[35 -15],backgroundImage);
hold on
% 
% if max(r_vec) > 8000
% brennend_image = importdata('brennendes_haus.jpg');
% %select the axes
% %axes(handles.Image);
% %place image onto the axes
% image([-15 -23 ],[3 -3],brennend_image)
% 
% else max(r_vec) <= 8000
%  nichtbrennend_image = importdata('nicht_brennend1.jpg');
% %select the axes
% %axes(handles.Image);
% %place image onto the axes
% image([-15 -23 ],[3 -3],nichtbrennend_image)
% 
% end
    


hold on

%-------------------- COlORBAR ----------------------

% COLORBAR: Image zu Beginn kleiner machen; Colorbar: Johannes und Jonas


min=100;
max=1200;
cbh=colorbar('location','eastoutside');
set(gca, 'CLim',[min,max]);
set(cbh,'XTickLabel',{num2str(min),num2str(max)})
% xlabel(cbh,'°C')



 
  z_vec = (hstep:hstep:z)';
 


 r_vec=r_vec/1000;
 z_vec=z_vec/1000;
 t_vec=t_vec/1000;
 
  plot (r_vec(r_vec~=0),z_vec(r_vec~=0),'r-')
  plot (-r_vec(r_vec~=0),z_vec(r_vec~=0),'r-')
 % pause(0.5)
  z2 =[z_vec; z_vec(end:-1:1)];
  r2= [r_vec;-r_vec(end:-1:1)];
 
 %t3=[t_vec;t_vec(end:-1:1)];
 t3=[t_vec(end:-1:1);t_vec];

  hold all
 % axis ([-35000 35000 -5000 35000]);
  axis([-40000/1000 40000/1000 -15000/1000 35000/1000])
  %set(gca,'xTickLabel',{'30','20','10','0','10','20','30'} )
  hold on
  xlabel('Breite (km)')
  ylabel('Höhe über Kraterradius (km)')
  set(gca,'xtick',[-40,-30,-20,-10,0,10,20,30,40],'xLim',[-40 40],...
      'XTickLabel',{'40','30','20','10',0,'10','20','30','40'});
  
  % ------------   Der Elipsenabschluss beim Plume  -----------------
  
  indx=find(r_vec);  % bringt mir nur die positiven Werte von r_vec hervor
  z_indx = z_vec(indx);
  z_plot = z_indx(end);
  t=linspace(0,2*pi,500); %von 0 bis 2*pi mit 500 schritten
  x= 0;
  y=z_plot;
  a=r_vec;
  b=1/10*r_vec; 
  plot(x-a*cos(t),y-b*sin(t),'Color',[.65 0 0]);  %dunkelrot
    indxnull= r_vec==0;  % sucht mir alle werte für r_vec=0 raus
    
%plot(z_vec(indxnull),t_vec(indxnull),'color','white' )
%t4=t3(r_vec~=0);
r3=r2(r2>0);
z4=z2(r2~=0);
r4=r2(r2~=0);


 hold on
 
%  -------------  Färbung des Kegels ------------------------------

 patch(r2(r2~=0),z2(r2~=0),t3(r2~=0),'edgecolor','none');

 hold on
% plot(z_vec(indxnull),'color','white' );
% z_vec(r_vec==0)=0;
%  pl=plot(z_vec(indxnull))
%  set(pl,'color','white')
%plot(z2(indxnull),t3(indxnull),'color','white' );% ich habe eigentlich einen langen...
  % strich bei r=0, von z=-50000 bis z_ende. Um den wegzubekommen plotte
  % ich alle z_vec-werte, die bei r_vec=0 sind weiss. so ist der strich
  % weg.
% 
 
 
 set(gca,'YDir','normal')  % beim Anzeigen im Gui wird das Bild umgedreht;
 % um das zu verhindern muss die Y-Achse hier angepasst werden fürs
 % GUI-Bild
   
 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
    if m11<=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));
    w1=(a1*y(4)*c1);
    if w1 <=0
        dy(1)=0;
    end
    dy(2)=2*al*rho_av*phi_av*(sqrt(a1 *y(4)*c1)) - omega*y(2)* d1;
     if w1 <=0
        dy(2)=0;
    end
    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)); 
        if w1 <=0
        dy(5)=0;
    end
    
   
   
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