%function diff_step
% solves the 1D diffusion problem for a step surface forcing
% animates the result
%
% Default parameter values at start are suitable for daily water:
%    Parameter values may be changed interactively
%    New values are saved for subsequent runs.
%    commands: "clear global" or "clear all" restore defaults
%

function diff_step
warning off MATLAB:divideByZero

global deltaT Tzero conductivity Cp rho Zmax tmax

clf

%default values (Water):
   if isempty(deltaT)
      deltaT=1;
   end
   if isempty(Tzero)
      Tzero=0;
   end
   if isempty(conductivity)
      conductivity=.600;
   end
   if isempty(rho)
      rho=1000;
   end
   if isempty(Cp)
      Cp=4200;
   end
   if isempty(Zmax)
      Zmax=1;
   end
   if isempty(tmax)
      tmax=1e4;
   end

change_param=menu('Change parameters?','no','chose material','enter values');
if change_param==2
   type=menu('Materials','liquid water','ice','air','rock','aluminum');
   if type==1
     %water
     conductivity= .603;
     Cp= 4200;
     rho=1000;
   elseif type==2
     %ice
     conductivity= 2.2;
     Cp= 2100;
     rho=1000;
   elseif type==3
     %air
     conductivity= 2.4e-2;
     Cp= 1000;
     rho=1.2;
   elseif type==4
     %rock
     conductivity= 2;
     Cp= 800;
     rho=2700;
   elseif type==5
     %aluminun
     conductivity= 233;
     Cp= 963;
     rho=2700;
   end
elseif change_param==3
  done = 0;
  while ~done
    k=menu('Choose parameter to enter','Delta T','Tzero',...
            'thermal conductivity','density','heat capacity',...
            'max depth','max time','none');
    if k==1
       deltaT=input(['Delta T = ',num2str(deltaT),...
                            ' Enter new value (^oC): ']);
    elseif k==2
       Tzero=input(['Tzero = ',num2str(Tzero),' Enter new value (^oC): ']);
    elseif k==3
       conductivity=input(['conductivity = ',num2str(conductivity),...
                            ' Enter new value (W/m/^oC): ']);
    elseif k==4
       rho=input(['density = ',num2str(rho),' Enter new value (kg/m^3: ']);
    elseif k==5
       Cp=input(['heat cpacity = ',num2str(Cp),...
                ' Enter new value (J/kg/^oC): ']);
    elseif k==6
       Zmax=input(['max depth = ',num2str(Zmax),' Enter new value (m): ']);
    elseif k==7
       tmax=input(['tmax = ',num2str(tmax),' Enter new value (s): ']);
    else
      done=1;
    end
  end
end


kappa = conductivity/(rho*Cp);
t=0:tmax/100:tmax;
z=0:Zmax/100:Zmax;

H211=subplot('position',[0.1300    0.7611    0.7750    0.1639]);
  cla
  set(gca,'ytick',[])
  set(gca,'xtick',[])
  set(gca,'box','on')
  text(.05,.85,'Current Parameters','fontweight','bold')

  text(.1,.67,['thermal conductivity = ',num2str(conductivity),' W/m/^oC'])
  text(.1,.48,['heat capacity = ',num2str(Cp),' J/kg/^oC'])
  text(.1,.29,['density = ',num2str(rho),' kg/m^3'])
  text(.1,.1,['diffusivity = ',num2str(kappa),' m^2/s'])

  text(.7,.67,['deltaT = ',num2str(deltaT),' ^oC'])
  text(.7,.48,['Tzero = ',num2str(Tzero),' ^oC'])
  text(.7,.29,['max depth = ',num2str(Zmax),' m'])
  text(.7,.1,['tmax = ',num2str(tmax),' s'])


disp(' ')
disp('Parameters')
disp(['> deltaT = ',num2str(deltaT),' degC'])
disp(['> Tzero = ',num2str(Tzero),' degC'])
disp(['> thermal conductivity = ',num2str(conductivity),' W/m/degC'])
disp(['> heat capacity = ',num2str(Cp),' J/kg/degC'])
disp(['> density = ',num2str(rho),' kg/m^3'])
disp(['> diffusivity = ',num2str(kappa),' m^2/s'])
disp(['> max depth = ',num2str(Zmax),' m'])
disp(['> tmax = ',num2str(tmax),' s'])
disp(' ')

Nt=max(size(t));
Nz=max(size(z));

Temperature=zeros(Nt,Nz);

for k=1:Nt
  Temperature(k,:)=Tzero+deltaT*erfc(z/(2*sqrt(kappa*t(k))) );
end

%subplot(121)
H121=subplot('position',[0.1300    0.1100    0.3270    0.6150]);
for k=1:Nt
  plot(Temperature(k,:),z)
  set(gca,'Ydir','reverse')
  xlabel('Temperature')
  ylabel('Depth (m)')
  text(Tzero+.5*deltaT,.8*Zmax,['time = ',num2str(t(k)),' s'])
  axis([Tzero Tzero+deltaT 0 Zmax])
  pause(.01)
end

done=0;
Tt_exists=0;
time=[];
Tt=[];

while ~done
  choice=menu('Choose depth to plot T(t)?','yes','no');
  if choice==2
    hold off
    done=1;
  else

  subplot(H121)
  %disp('Pick depth with the mouse')
  text(Tzero+.1*deltaT,.9*Zmax,['Pick depth with the mouse'],'fontsize',12)
  Tt_exists=1;
  
  [dummy,depth]=ginput(1);
    if depth < 0
       depth=0;
    end
    if depth > Zmax
       depth=Zmax;
    end
  %Nz
  %depth
  %Zmax
  N=round(Nz*depth/Zmax);
  text(Tzero+.1*deltaT,.9*Zmax,['Pick depth with the mouse'],...
         'fontsize',12,'color','w')

  H122=subplot('position',[0.5780    0.1100    0.3270    0.6150]);
    plot(t,Temperature(:,N))
    ylabel('Temperature (^oC)')
    xlabel('Time (s)')
    axis([0 tmax Tzero Tzero+deltaT])
    text(.5*tmax,Temperature(round(.5*Nt),N),['depth = ',num2str(z(N)),' m'])
    hold on

  end
end

if Tt_exists
done=0;

while ~done
  choice=menu('Pick a point on T(t) curve?','yes','no');
  if choice==1
    subplot(H122)
    if ~isempty(time) & ~isempty(Tt)
      text(0.1*tmax,Tzero+0.9*deltaT,...
            ['t = ',num2str(time,3),'s'],'color','w')
      text(0.1*tmax,Tzero+0.8*deltaT,...
            ['T = ',num2str(Tt,3),' ^oC'],'color','w')
    end
    text(0.1*tmax,Tzero+0.9*deltaT,'Pick point with the mouse','fontsize',12)
    point=pickxy;
    text(0.1*tmax,Tzero+0.9*deltaT,...
          'Pick point with the mouse','fontsize',12,'color','w')
    time=point(1);
    Tt=point(2);
      hold on
      ylabel('Temperature (^oC)')
      xlabel('Time (s)')
      axis([0 tmax Tzero Tzero+deltaT])
      text(.5*tmax,Temperature(round(.5*Nt),N),['depth = ',num2str(z(N)),' m'])
      text(0.1*tmax,Tzero+0.9*deltaT,['t = ',num2str(time,3),'s'])
      text(0.1*tmax,Tzero+0.8*deltaT,['T = ',num2str(Tt,3),' ^oC'])
      plot(time,Tt,'go','markerfacecolor','g')
  else
    done=1;
  end
end

end



