% Gorelick for GES 231/CEE260C      April 9,2001/02/03
% Advection Dispersion Eqn in 1-D  : Analytic1AD_movie.m
% Analytic Solution -- Semi-Infinite Case
% Fixed Concentration BC on Left; Zero at infinity on right
% Modified by Trayle Kulshan April 10, 2002: added movie

% Modified by Nick Martin April 9, 2003: fix divide by zero and add
% movie switch; MovieOn=1 plays movie and MovieOn=0 for not movie.

% NOTE: This program requires that experfc.m be in the same directory

% Produces Concentrations Profiles at 3 specified times (FIG 2)
% Produces Concentration Histories at 3 specified locations (FIG 3)
% User must adjust labels (hard wired) on Figure 2 and 3 as needed
% Self consistent units are used for all constants and variables
% Concentrations are in arbitrary units (system linear in concentration)

% identify constants
% clear lines modified for "new" syntax. NM 04/09/2003
clear all;
close all hidden;
fidprofiles= fopen('analytic1prof.out','w+');
fidhistories= fopen('analytic1hist.out','w+');

%NOTE -- SET q>0 below to generate solution for ADVECTIVE-Dispersive Transport
%you will have to do this for Homework Assignment 

MovieOn = 1;    % 1 means that movie will play; 0 means that movie will not play.
por=0.3;        % effective porosity
q=0.0;          % specific discharge, m/d <- MODIFY THIS LINE TO INCLUDE ADVECTION
DH=2.5;         % hydrodynamic disp. coef = alfa * |q|= alfa*por*|V|
V=q/por;        % groundwater velocity --> average linear gw velocity
DHP=DH/por;     % hydrodynamic disp. coef. = alfa * V
%NOTE: constants V and DHP used here do not include effective porosity; 
% if your numerical model is based on an equation that contains 
% effective porosity on left,then divide through to get V and DHP used here 
% if you want to compare to this solution.

cbound=100.0; % concentration at boundary on left

% space and time paramters for concentration profiles
length=50;       % Length over which results are produced 
DX=1.0;     % distance increment at which results are produced 
N=(length/DX)+1; % number of nodes used in numerical model at which
                                     % analytic solution is to be computed
node_locations=[0:DX:length]; % vector containing position of nodes
times=[1.0 3.0 5.0];          % times at which concentration profiles are produced

% time and space markers for concentration histories
time_frame=5.0;  % time horizon of problem (days)
DT=0.05;         % time increment
number_timesteps=(time_frame/DT)+1; % number of time steps in numerical model at 
                                    % which analytic solution is to be computed

fixed_locations=[4.0 14.0 20.0];  % locations at which conc histories are 
                                  % produced;DISTANCES from left boundary 
time_increments=[0:DT:time_frame];% vector containing all time steps

% Begin production of concentration profile movie
% TIME LOOP produced profiles at every timestep
for t=DT:DT:time_frame
   % SPACE LOOP produces profiles at nodal locations,ii
    c(1)=cbound;
    for ii=2:1:N;       
        x=node_locations(ii);
        % lines added by NM 04/09/2003 to fix divide by zero.%%%%
        temp1 = DHP;
        btemp = (abs(temp1 - double(0.0)) > eps);
        temp2 = (btemp*temp1) + ((1 - btemp)*eps);
        temp1 = 1/temp2;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        arg1=(V*x)*temp1;
        % lines added by NM 04/09/2003 to fix divide by zero.%%%%
        temp1 = (2*(sqrt(DHP*t)));
        btemp = (abs(temp1 - double(0.0)) > eps);
        temp2 = (btemp*temp1) + ((1 - btemp)*eps);
        temp1 = 1/temp2;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        arg2=(x+V*t)*temp1;
        arg3=experfc(arg1,arg2);
        c(ii)=(cbound/2)*((erfc((x-V*t)*temp1))+arg3);
     end
     
%Swith to turn of movie is checked here. Use set value above.     
    if (MovieOn == 1)
        figure(1);
        clf;
        hold on;
        plot(node_locations,c,'b-o');
        xlabel('distance [m]');
        ylabel('Concentration [mg/L]');
        title(['Concentration Profile - Analytic Soln "O" at time: ' num2str(t,'%4.2f') ' days']);
        axis([0 50 0 100]);
        pause(0.01);
        hold off;
    end
end
% End of movie 

% Begin production of concentration profiles at specified times
% TIME LOOP produced profiles at 3 times, t
for timeindex=1:1:3;  
   t=times(timeindex);
   % SPACE LOOP produces profiles at nodal locations,ii
   c(1)=cbound;
   for ii=2:1:N;        
      x=node_locations(ii);
      % lines added by NM 04/09/2003 to fix divide by zero.%%%%
      temp1 = DHP;
      btemp = (abs(temp1 - double(0.0)) > eps);
      temp2 = (btemp*temp1) + ((1 - btemp)*eps);
      temp1 = 1/temp2;
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      arg1=(V*x)*temp1;
       % lines added by NM 04/09/2003 to fix divide by zero.%%%%
      temp1 = (2*(sqrt(DHP*t)));
      btemp = (abs(temp1 - double(0.0)) > eps);
      temp2 = (btemp*temp1) + ((1 - btemp)*eps);
      temp1 = 1/temp2;
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      arg2=(x+V*t)*temp1;
      arg3=experfc(arg1,arg2);
      c(ii)=(cbound/2)*((erfc((x-V*t)*temp1))+arg3);
   end
      %store 3 concentration profiles at specified times (days 1,3,5)

      if t==times(1)
        c1=c;
             c1;
          elseif t==times(2) 
          c3=c;

          c3;
          elseif t==times(3)  
          c5=c;
          c5;
        end
end
     
% store 3 concentration profiles c(x)=locations,concentrations 

fprintf(fidprofiles,'%.6f\t',node_locations);
fprintf(fidprofiles,'\n');

disp ('Locations and Concentrations for Time=1');
c1_loc_and_conc=[node_locations; c1]
fprintf(fidprofiles,'%.6f\t',c1);
fprintf(fidprofiles,'\n');

disp ('Locations and Concentrations for Time=3');
c3_loc_and_conc=[node_locations; c3]
fprintf(fidprofiles,'%.6f\t',c3);
fprintf(fidprofiles,'\n');

disp ('Locations and Concentrations for Time=5');
c5_loc_and_conc=[node_locations; c5]
fprintf(fidprofiles,'%.6f\t',c5);
fprintf(fidprofiles,'\n');
fclose(fidprofiles);
% plot figure with concentration profiles -- label locations are fixed
% by user using crosshairs and mouse click
figure(2);
hold on;
plot(node_locations,c1,'-+',node_locations,c3,'-+',node_locations,c5,'-+');
xlabel('distance [m]');

ylabel('Concentration [mg/L]');
title('Concentration Profiles - Analytic Soln "+"');
gtext('Day 1');gtext('Day 3');gtext('Day 5');
hold off;

%code below tests overwriting of analytic solution on Figure 1
%  hold on; 

%  load('analytic1prof.out');
%  J2=analytic1prof(1,:);
%  plot(J2,analytic1prof(2:4,:),'o');
% End of concentration profiles at specified times

% Begin production of concentration histories at fixed locations
%  Space LOOP produced histories at 3 locations, x
for location_index=1:1:3;  
   x=fixed_locations(location_index);
   % time LOOP produces histories at times jj 
   for jj=1:1:number_timesteps; 
      tt=time_increments(jj);

      % lines added by NM 04/09/2003 to fix divide by zero.%%%%
      temp1 = DHP;
      btemp = (abs(temp1 - double(0.0)) > eps);
      temp2 = (btemp*temp1) + ((1 - btemp)*eps);
      temp1 = 1/temp2;
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      arg1=(V*x)*temp1;
      % lines added by NM 04/09/2003 to fix divide by zero.%%%%
      temp1 = 2*(sqrt(DHP*tt));
      btemp = (abs(temp1 - double(0.0)) > eps);
      temp2 = (btemp*temp1) + ((1 - btemp)*eps);
      temp1 = 1/temp2;
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      arg2=(x+V*tt)*temp1;
      arg3=experfc(arg1,arg2);
      cc(jj)=(cbound/2)*((erfc((x-V*tt)*temp1))+arg3);
   end
   %store 3 concentration histories at specified locations
        cc(1)=0.0;
        if x==fixed_locations(1);
        c1_hist=cc;
             c1_hist;
          elseif x==fixed_locations(2);   
          c2_hist=cc;
        c2_hist;
        elseif x==fixed_locations(3); 
          c3_hist=cc;
          c3_hist;
                end
end

% store 3 concentration histories c(t)=locations,concentrations 

disp ('Times and Concentrations for Location 1');
c1_time_and_conc=[time_increments; c1_hist]
fprintf(fidhistories,'%.6f\t',time_increments);
fprintf(fidhistories,'\n');
fprintf(fidhistories,'%.6f\t',c1_hist);
fprintf(fidhistories,'\n');


disp ('Times and Concentrations for Location 2');
c2_time_and_conc=[time_increments; c2_hist]
fprintf(fidhistories,'%.6f\t',c2_hist);

fprintf(fidhistories,'\n');


disp ('Times and Concentrations for Location 3');
c3_time_and_conc=[time_increments; c3_hist]
fprintf(fidhistories,'%.6f\t',c3_hist);
fprintf(fidhistories,'\n');

% plot figure with concentration profiles -- label locations are fixed
% by user using crosshairs and mouse click
figure(3);
hold on;
plot(time_increments,c1_hist,'-+',time_increments,c2_hist,'-+',...
   time_increments,c3_hist,'-+');
xlabel('time [days]');
ylabel('Concentration [mg/L]');
title('Concentration Histories - Analytic Soln "+"');
gtext('Location 4 m');gtext('Location 14 m');gtext('Location 20 m');
hold off;

fclose(fidhistories);

%code below tests overwriting of analytic solution on Figure 2
%  hold on; 
%  load('forwarddif_histories.out');
%  J2=analytic1hist(1,:);
%  plot(J2,analytic1hist(2:4,:),'o');
% End of concetration histories at fixed locations

fclose('all');

%EOF
