% The model describes a scenario of multiple intrusions (in 1950, 1970,
% 1982 and continuous intrusions since 2005) from 1950 until 2035 in a sill
% of 25 m or 15 m in thickness for two volumes estimated from geodetic
% inversion

clc
clear all

mfile_name          = mfilename('fullpath');
[pathstr,name,ext]  = fileparts(mfile_name);
cd(pathstr);

fileposition_model    = 'model_outputs/25m_volmax/' %change the last 
%folder's name to explore the other models (25m_volmin/, 25,_volmax/, 
%15_volmax/ or 15m_volmin/). Do not delete "model_outputs/".

load([fileposition_model,'setup']) 

it=1:length(time_yrs)
   
fraction(it) = length(find(displ(kii)>=Teruptible))/length(kii);


for n = 850 
xmax=5;

   load([fileposition_model,'output_',num2str(n)])  
   load([fileposition_model,'output_t']) 

   solution   = reshape(displ,nx,ny, [])' ;

   figure(1)  ; clf
   pcolor(xgrid,ygrid,solution)  
   colorbar
   xlabel('Horizontal distance (m)')
   ylabel('Depth (m)')
   title(['Intrusion after ',num2str(time(n)/seconds_per_year),' years'])
   shading interp

   figure(2)
   plot(time/seconds_per_year,meanT,'b')
   xlabel('Time (y) from 1950')
   ylabel('Temperature (°C)')
   title('Mean sill temperature (°C) vs time (y)')
   xlim([0 85])

   figure(3) %Histogram of intrusion temperatures
   hist(particles_T_t(:,it))
   title('Histogram of intrusion temperatures (°C)')


   figure(4) %Temperature profile through the center of the sill
   plot(displ(bx0),g_coord(2,bx0)/1e3)   
   xlabel('Temperature (°C)')
   ylabel('Depth (km)')
   title('Temperature vs depth (km) through center of intrusion')
   axis([0 Tintrusion -ly/1e3 -0])


   figure(5) %area of magma for each isotherms versus time
   plot(time/seconds_per_year,areas*1e-6)
   title('Area greater than each isoterm (km^2) vs time (y)')
   ylabel('Area for each isotherm (km^2)')
   xlabel('Time (y) from 1950')
   title(['Area of magma for a sill of h=',num2str(Ithickness),' meters'])
   leg = legend('show'); 
   title(leg,'Temperature (°C)')
   legend('>850', '>900', '>950','>1000', '>1050', '>1100','>1150','>1200','>1250')
   xlim([0 85])
  
   volume_km3 = vol*1e-9;

   figure(6) %volume of magma for each isotherms versus time
   plot(time/seconds_per_year,volume_km3)
   title('Volume greater than each isoterm (km^3) vs time (y)')
   ylabel('Volume for each isotherm (km^3)')
   xlabel('Time (y) from 1950')
   title(['Volume of magma for a sill of h=',num2str(Ithickness),' meters'])
   leg = legend('show'); 
   title(leg,'Temperature (°C)')
   legend('>850', '>900', '>950','>1000', '>1050', '>1100','>1150','>1200','>1250')
   xlim([0 85])



end