%function diff_sin
% solves the 1D diffusion problem for a sinusoidal surface forcing
% animates the result
%
% Default parameter values at start are suitable for daily variation and water:
%    Parameter values may be changed interactively
%    New values are saved for subsequent runs.
%    commands: "clear global" or "clear all" restore defaults
%
function diff_sin
warning off MATLAB:divideByZero
                                                                                
global deltaT Tzero conductivity Cp rho Zmax period
                                                                                
%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=.3;
   end
   if isempty(period)
      period=86400;
   end

%while working
  change_param=menu('Change parameters?','no','chose material','enter values');
  if change_param==2
     type=menu('Material:','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','Peak-to-peak delta T','Tzero',...
              'thermal conductivity','density','heat capacity',...
              'max depth','period','none');
      if k==1
         deltaT=input(['Delta T = ',num2str(deltaT),...
                              ' Enter new value (deg C): ']);
      elseif k==2
         Tzero=input(['Tzero = ',num2str(Tzero),' Enter new value (deg C): ']);
      elseif k==3
         conductivity=input(['conductivity = ',num2str(conductivity),...
                              ' Enter new value (W/m/degC): ']);
      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/deg C): ']);
      elseif k==6
         Zmax=input(['max depth = ',num2str(Zmax),' Enter new value (m): ']);
      elseif k==7
         period=input(['period = ',num2str(period),' Enter new value (s): ']);
      else
        done=1;
      end
    end
  end
                                                                                
  omega=2*pi/period;
  kappa = conductivity/(rho*Cp);
  delta=sqrt(2*kappa/omega);
  
  subplot(2,2,4)
    cla
    set(gca,'ytick',[])
    set(gca,'xtick',[])
    set(gca,'box','on')
  disp(' ')
  disp('Parameters')
    text(.025,.9,'Current Parameters','fontweight','bold')
  
  disp(['> peak-to-peak deltaT = ',num2str(deltaT),' degC'])
    text(.04,.8,['peak-to-peak deltaT = ',num2str(deltaT),' degC'])
  
  disp(['> Tzero = ',num2str(Tzero),' degC'])
    text(.04,.7,['Tzero = ',num2str(Tzero),' degC'])
  
  disp(['> thermal conductivity = ',num2str(conductivity),' W/m/degC'])
    text(.04,.6,['thermal conductivity = ',num2str(conductivity),' W/m/degC'])
  
  disp(['> heat capacity = ',num2str(Cp),' J/kg/degC'])
    text(.04,.5,['heat capacity = ',num2str(Cp),' J/kg/degC'])
  
  disp(['> density = ',num2str(rho),' kg/m^3'])
    text(.04,.4,['density = ',num2str(rho),' kg/m^3'])
  
  disp(['> diffusivity = ',num2str(kappa),' m^2/s'])
    text(.04,.3,['diffusivity = ',num2str(kappa),' m^2/s'])
  
  disp(['> max depth = ',num2str(Zmax),' m'])
    text(.04,.2,['max depth = ',num2str(Zmax),' m'])
  
  disp(['> period = ',num2str(period),' s'])
    text(.04,.1,['period = ',num2str(period),' s'])
  
  disp(' ')
  
  rad=180/pi;
  while(1)
    z=0:Zmax/100:Zmax;
    tmax=2*period;
    N=100;
    t=0:tmax/N:tmax;
  
    temp0= Tzero+0.5*deltaT*imag(exp(i*omega*t));
    amp=deltaT*exp(-z/delta);

    subplot(2,2,3)
      hold off
      plot(amp,z,'b','linewidth',2)
      axis([0 deltaT 0 Zmax])
      ylabel('Depth (m)')
      xlabel('Temperature variation (degC)')
      set(gca,'Ydir','reverse')
      kk=menu('Change maximum depth?','yes','no')
      if kk==1
        Zmax=input('Input Zmax (m): ')
      else
        break
      end
    end
    
  
  subplot(2,1,1)
  hold off
    for k=1:N
      %plot(t(k),temp0(k),'ro')
      plot(t(1:k),temp0(1:k),'r','linewidth',2)
      hold on
      axis([0 tmax Tzero-deltaT/2 Tzero+deltaT/2])
      title('Time variation of surface temperature')
      ylabel('Temperature (degC')
      xlabel('Time (s)')
    pause(1e-10)
  end
  
  exist_lagged=0;
  done=0;
                                                                                  
while ~done
    choice=menu('Choose depth to plot T(t)?','yes','no');
    if choice==2
      hold off
      done=1;
    else
      while(1)
        exist_lagged=1;
          subplot(2,2,3)
          hold off
          z=0:Zmax/100:Zmax;
          amp=deltaT*exp(-z/delta);
          plot(amp,z,'linewidth',2)
          axis([0 deltaT 0 Zmax])
          ylabel('Depth (m)')
          xlabel('Temperature variation (degC)')
          set(gca,'Ydir','reverse')
          text(0.1*deltaT,0.9*Zmax,'Pick a depth with the mouse','fontsize',12)
          subplot(2,2,3)
          [dummy depth]=ginput(1);
            if depth < 0
               depth=0;
            end
            if depth > Zmax
               depth=Zmax;
            end
  
          ampZ=deltaT*exp(-depth/delta);
          tempZ= Tzero+0.5*ampZ*imag(exp( i*(omega*t - depth/delta) ));
          hold off
          plot(amp,z,'b','linewidth',2)
          axis([0 deltaT 0 Zmax])
          ylabel('Depth (m)')
          xlabel('Temperature variation (degC)')
          set(gca,'Ydir','reverse')
          hold on
          plot(ampZ,depth,'gd','markerfacecolor','g')
%          plot(ampZ,depth,'g','linewidth',2)
          text(0.58*deltaT,0.8*Zmax,['depth = ',num2str(depth,2)])
          text(0.58*deltaT,0.9*Zmax,['amplitude = ',num2str(ampZ,2)])
          kk=menu('Change maximum depth?','yes','no')
          if kk==1
            Zmax=input('Input Zmax (m): ')
          else
            break
          end
      end
  
      subplot(2,1,1)
        hold off
        for k=1:N
 %         plot(t(k),temp0(k),'ro')
          plot(t(1:k),temp0(1:k),'r','linewidth',2)
          hold on
%          plot(t(k),tempZ(k),'gd')
          plot(t(1:k),tempZ(1:k),'g','linewidth',2)
          axis([0 tmax Tzero-deltaT/2 Tzero+deltaT/2])
          title('Time variation of temperature')
          ylabel('Temperature (degC)')
          xlabel('Time (s)')
          pause(1e-10)
        end
        plot([0 tmax],[Tzero Tzero])
        plot([tmax/2 tmax/2],[Tzero-deltaT/2 Tzero+deltaT/2])
    end
end
  
if exist_lagged
  done=0;
  while ~done
    choice=menu('Pick zero-crossing?','yes','no');
    if choice==2
      hold off
      done=1;
    else
      subplot(2,1,1)
        hold off
        %plot(t,temp0,'ro')
        plot(t,temp0,'r','linewidth',2)
        hold on
        %plot(t,tempZ,'gd')
        plot(t,tempZ,'g','linewidth',2)
        axis([0 tmax Tzero-deltaT/2 Tzero+deltaT/2])
        title('Time variation of temperature')
        ylabel('Temperature (degC)')
        xlabel('Time (s)')
        plot([0 tmax],[Tzero Tzero])
        plot([tmax/2 tmax/2],[Tzero-deltaT/2 Tzero+deltaT/2])
        time=pickx;
        plot([time time],[Tzero-deltaT/2 Tzero+deltaT/2],'g--')
        hold off
        text(.75*tmax,(Tzero+.8*deltaT/2),...
             ['phase lag = ',num2str(time-tmax/2,3)]) 
    end
  end
end
