% Gorelick Example Program             April 10, 2001
% Forward Difference Model for
% Dispersion Eqn in 1-D
% Basis for Homework 1 (a) Run as is and then ADD Advective Term

% Open files for concentration profiles and concentration histories
clear variables

% identify constants

DT=.05;     % Time step
DX=1;		   % Node interval
por=0.3;		% Effective Porosity
Cf=100.0;	% Influx Conc. at x=0 t>0  [mg/L]
Cb=0.0;		% Boundary cond at x=50 m and initial cond at all x
DH=2.5;		% hydrodynamic dispersion coef. [m2/d]= = alfa * |q|= alfa*por*|V|

L=50;			% Length of system [m]
Tmax=5;		% Max time of simulation
N=L/DX+1;   % Number of Nodes

node_locations=[0:DX:L]; % vector containing position of N nodes
J=[0:DX:L];	% shorthand for vector containing position of N nodes
fixed_locations=[5.0 15.0 21.0];  % locations at which conc histories are 
											 % produced;DISTANCES from left boundary 
number_timesteps=(Tmax/DT); % number of time steps in numerical model at 
									 % which analytic solution is to be computed
time_increments=[0:DT:Tmax];% vector containing all time steps

%Initialize Concentrations and Boundary Values
Cold=ones(1,N)*Cb;	% C is a row vector; Initialize for t=0
Cold(1)=Cf;
Cnew=Cold;
% Set coefficients in FD equations
p=(DH*DT)/(por*DX^2); 	

a=p;		
b=1.0-2*p;		
c=p;		

% Time step loop--marching through time at DT time steps

NPRINT=0;

% 	Open output file and write the initial parameters
%	Boundary, Initial Conditions; Time step; Node Interval;
%	Dispersion Coefficient; GW flux
fid= fopen('fordif.out','w+');
fprintf(fid,'\tLeft Boundary Condition Cf=%.3f [mg/L]\n',Cf);
fprintf(fid,'\tRight Boundary Condition Cb=%.3f [mg/L]\n',Cb);
fprintf(fid,'\tInitial Condition at x=0  : C(0)=%.3f [mg/L]\n',Cf);
fprintf(fid,'\tInitial Condition at x>0  : C(i)=%.3f [mg/L]\n',Cb);
%fprintf(fid,'\tGw flux q=%.3f [m/d]\n ',q);
fprintf(fid,'\tDisp.Coef. DH=%.3f [m2/d]\n ',DH);
fprintf(fid,'\tTime Step DT=%.3f [days]\n',DT);
fprintf(fid,'\tNode Interval DX=%.3f [m]\n',DX);
fprintf(fid,'\tTotal Nodes N=%1g\n',N);
fprintf(fid,'\tSystem Length L=%.3f [m]\n',L);
fprintf(fid,'\tMax Time TMAX=%.3f [m]\n',Tmax);
fprintf(fid,'\n----------------------------\n Intermediate Parameters:');
fprintf(fid,'\np=\t%.4f',p);

fprintf(fid,'\na=\t%.4f',a);
fprintf(fid,'\nb=\t%.4f',b);
fprintf(fid,'\nc=\t%.4f\n',c);

% TIME LOOP by DT up to Tmax
% Space Loop goest from 1 to 50; note it does not start at zero
% and value at node 1 and 51 are known from boundary conditions

% Note how time step vector is st up
% time steps and time step vector
% 0  dt 2dt 3dt....tmax
% 1  2   3   4     time_increments
%Time set counter -- vector of time steps starts at itime=1 and is time zero
%Time step counter 2, vector location 2, is first time setp DT
itime=1;
for i=DT:DT:Tmax;    % i in days
   Cnew(1)=Cf;  	% Boundary condition on left 
   Cnew(51)=Cb;
   NPRINT=NPRINT+1;
	itime=itime+1;	
	for k=2:1:N-1;	% Nodes 2 to N-1 are actively solved for
      Cnew(k)=(a*Cold(k+1))+((b)*Cold(k))+(c*Cold(k-1));
%store 3 concentration histories at specified locations
       
		  if k==fixed_locations(1);
        c1_hist(itime)=Cnew(k);
   	  elseif k==fixed_locations(2);   
   	  c2_hist(itime)=Cnew(k);
        elseif k==fixed_locations(3);   
        c3_hist(itime)=Cnew(k);
        end
	end
		Cold=Cnew; 
	if NPRINT==10;
		fprintf(fid,'Time step=\t%.3f days\n\n',i);
		fprintf(fid,'%.6f\t',Cnew);
		fprintf(fid,'\n\n');
		NPRINT=0;
	end
	
%	Collect Profiles at days 1,3,5	
	if i==1
      C1=Cnew;
	elseif i==3 
      C3=Cnew;
   elseif i==5
      C5=Cnew;
   end
end
% store 3 concentration profiles c(x)=locations,concentrations 
fidprofiles= fopen('fordif_profiles.out','w+');
fidhistories= fopen('fordif_histories.out','w+');

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);

% 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');

%	Forward Difference Solution noted with "+" Concentration Profile
figure(3);
plot(J,C1,'+',J,C3,'+',J,C5,'+');
hold on;
xlabel('distance [m]')
ylabel('Concentration [mg/L]')
title('HW1a- Forward Difference Soln vs Analytic Soln: Profiles')
legend('Numerical Soln')
gtext('Day 1');gtext('Day 3');gtext('Day 5');

%code below overwrites analytic solution on Figure 1 - Conc. Profiles

hold on; 
load('analytic1prof.out');
J2=analytic1prof(1,:);
plot(J2,analytic1prof(2:4,:),'-o');

%  Forward difference solution Conc. histories
figure(4);
plot(time_increments,c1_hist,'+',time_increments,c2_hist,'+',...
   time_increments,c3_hist,'+');
xlabel('time [d]')
ylabel('Concentration [mg/L]')
title('HW1a- Forward Difference Soln vs Analytic Soln: Histories')
legend('Numerical Soln')
gtext('Location 1');gtext('Location 2');gtext('Location 3');

%code below overwrites analytic solution on Figure 2- Con. Histories
hold on;
load('analytic1hist.out');
J2=analytic1hist(1,:);
plot(J2,analytic1hist(2:4,:),'-o');

fclose('all');

 	



