function ComCZAR_RK4(parfilename)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  ComCZAR ---  A 2-D Homogeneous Anisotropic Complex Potential Analytical Model
%           for Capture Zone Delineation.  Arbitrary number of wells with
%           Arbitrary flow rates, and optional regional flow.
%       Created by Mike Fienen, Jian Luo, and Peter Kitanidis
%       Stanford University - Department of Civil and Environmental Engineering
%       Environmental Fluid Mechanics Laboratory
%   This code accompanies a paper in revision for Journal of Hydrology (12/2004) entitled:
%       "Semi-Analytical Homogeneous Anisotropic Capture Zone Delineation"
%          
%   This version of the code uses a 4th order Runge-Kutta scheme
%   	for streamline tracing.
%   The primary contact with bug reports is Mike Fienen (fienen@stanford.edu)
%   Also, if you use this code, reference the paper above, and please drop us a 
%       note letting us know how it worked.
%   Further documentation can be found at http://www.stanford.edu/~fienen/software/            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%READ IN PARAMETERS
    parfile = fopen(parfilename,'r');
	%skip header lines
    	for j = 1:6; junkus	= fgetl(parfile); end
    %get figure filename
        params_in = fgetl(parfile);
            params_in = sscanf(params_in,'%s');
            for j = 1:length(params_in)
                if params_in(j) ~= '*'
                    outfigname(j) = params_in(j);
                else break
                end
            end
         outfigname = ['RK4' outfigname];
         params_in	= fgetl(parfile);
            num_iter = sscanf(params_in, '%f');
    %read anisotropy parameters
        junkus	= fgetl(parfile); %skip explanation line
        params_in	= fgetl(parfile);
            Kmax_Kmin	= sscanf(params_in, '%f'); 
        params_in	= fgetl(parfile);
            theta	= sscanf(params_in, '%f');
    %read regional flow parameters
        junkus	= fgetl(parfile); %skip explanation line
        params_in	= fgetl(parfile);
            T   	= sscanf(params_in, '%f'); 
        params_in	= fgetl(parfile);
            J   	= sscanf(params_in, '%f');
        params_in	= fgetl(parfile);
            gamma  	= sscanf(params_in, '%f');    
    %read discretization parameters
        junkus	= fgetl(parfile); %skip explanation line
        params_in	= fgetl(parfile);
            dx   	= sscanf(params_in, '%f');
        params_in	= fgetl(parfile);
            dy   	= sscanf(params_in, '%f');
        params_in	= fgetl(parfile);
            xmin   	= sscanf(params_in, '%f');
        params_in	= fgetl(parfile);
            xmax   	= sscanf(params_in, '%f');
        params_in	= fgetl(parfile);
            ymin   	= sscanf(params_in, '%f');
        params_in	= fgetl(parfile);
            ymax   	= sscanf(params_in, '%f');
    %skip explanation lines
    	for j = 1:6; junkus	= fgetl(parfile); end 
    %read well data
        j = 1;
        while (~feof(parfile))
            params_in	= fgetl(parfile);
            [wellname{j},wells(j,1),wells(j,2),Q(j)]=strread(params_in,'%s%f%f%f');
            j = j+1;
        end
        fclose(parfile);
%DONE WITH PARAMTERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    i       = sqrt(-1);

%Set the stretching coefficient and the angle of principle anisotropy
	bet     = sqrt(Kmax_Kmin);
	theta   = theta*pi/180;
        
%Create the rotation matrix
    RM      = [cos(theta) sin(theta);-sin(theta) cos(theta)];
%Convert the flow rates, initially in liters per minute, into m^3/sec
	Q       = Q./1000./60;

%discretize the domain for plotting potential
    x       = [xmin:dx:xmax];
	y       = [ymin:dy:ymax];
	[X,Y]   = meshgrid(x,y);
    Z       = X + i*Y;
    
%Transform the calculation domain for the analytical solution 
    X_tr    = X*cos(theta) + Y*sin(theta);
    Y_tr    = -X*sin(theta) + Y*cos(theta);

%stretch transform of the domain
    Z_tr    = X_tr + i*Y_tr*bet;
    
%Set k and determine regional flow.  If k==0, no regional flow.  This k value does
%not affect the rest of the solution.
    gamma   = gamma*pi/180;
    Tx       = T; % m^2/s
    Ty       = Tx/bet^2;
    Jz       = J*exp(i*gamma);
    Jz_t_tmp = RM*[real(Jz); imag(Jz)];
    Jz_t     = Jz_t_tmp(1) + i*Jz_t_tmp(2);
    Qr       = Tx*real(Jz_t)/bet + i*Ty*imag(Jz_t);  
%Transform the well locations
[lwells,junk] = size(wells);
    for j = 1:lwells
        wells_tr(j,:) = transpose(RM*wells(j,:)');
        wells_Z_tr(j) = wells_tr(j,1) + i*bet*wells_tr(j,2);
    end
    
%Calculate OMEGA, PHI, and PSI in transformed domain
    Omt_r = -conj(Qr)*Z_tr;
	for j = 1:lwells
        Omt_r = Omt_r + Q(j)/(2*pi)*log(Z_tr-wells_Z_tr(j)+eps);
	end
	phit_r = real(Omt_r);

    outfig = figure;
    hold on;
%Plot the PHI 
    [cs,ch]=contour(X,Y,phit_r, 50,'k--');

%free up some memory
    clear phi_tr Omt_r;

%Find stagnation points analytically    
    a = Q./2./pi;  ww=transpose(wells_Z_tr);
    ww(a==0) = []; a(a==0)=[];
    exw      = ww(a<0);
    [F, G]      = residue(-a, ww, [conj(Qr)]);
    z_stag_tr   = roots(F);

    stag_y_tr   = imag(z_stag_tr);
    stag_x_tr   = real(z_stag_tr);
    stag_x      = stag_x_tr*cos(-theta) + stag_y_tr/bet*sin(-theta);
    stag_y      = -stag_x_tr*sin(-theta) + stag_y_tr/bet*cos(-theta);

%Plot the stagnation streamlines   
    Qplot = Q;
    Q(Q==0) = [];
    zw_tr  =   transpose(ww);
    z_stag_part_maj = [];
    z_stag_part_min = [];
    ds              =   sqrt(dx^2 + dy^2);
    offset          =   ds/25;
%Contruct the Hessian for the flow field at each stagnation point
%Elemental formulae were derived analytically from the potential
for j = 1:length(z_stag_tr)

    pdz2_Om = sum(-Q./2./pi./(zw_tr-z_stag_tr(j) + eps).^2);      
    H=[real(pdz2_Om) -imag(pdz2_Om);-imag(pdz2_Om) -real(pdz2_Om)];
    
    
%Find the eigenvalues and eigenvectors of the Hessian matrix
    [e_vect,e_val] = eig(H);
    
    if e_val(1) > 0
        p_dir   = e_vect(:,1);    
        min_dir = e_vect(:,2);
    else 
        p_dir   = e_vect(:,2);
        min_dir = e_vect(:,1);
    end
%determine points in both the major and minor directions emanating from the stagnation points
%for particles to be tracked
    sx_maj = [stag_x_tr(j)+offset*p_dir(1); stag_x_tr(j)-offset*p_dir(1)];
    sy_maj = [stag_y_tr(j)+offset*p_dir(2); stag_y_tr(j)-offset*p_dir(2)];
    z_stag_part_maj = [z_stag_part_maj;sx_maj + i*sy_maj];    
    sx_min = [stag_x_tr(j)+offset*min_dir(1); stag_x_tr(j)-offset*min_dir(1)];
    sy_min = [stag_y_tr(j)+offset*min_dir(2); stag_y_tr(j)-offset*min_dir(2)];
    z_stag_part_min = [z_stag_part_min;sx_min + i*sy_min];   
end    
    pnum    =   length(z_stag_part_maj);

%TRACK THE PARTICLES 
    zp_maj      =   reshape(repmat([stag_x_tr + i*stag_y_tr].',2,1),length(stag_x_tr)*2,1);
    zp_min      =   zp_maj;
    zp_maj(:,1) =   z_stag_part_maj;
    zp_min(:,1) =   z_stag_part_min;

    dr_maj = min(abs(repmat(zp_maj(:,1).',length(ww),1) - repmat(ww,1,length(zp_maj))));
    dr_min = min(abs(repmat(zp_min(:,1).',length(ww),1) - repmat(ww,1,length(zp_min))));
	N = 2; 
    while ( N<=num_iter )
        [vz_maj,a,zp4_maj] = tracer_RK4(zp_maj(:,N-1), -ds, theta, bet, Q, ww, Qr);
        zp_maj(:,N)        = zp4_maj + ds.*exp(i*(a+pi));
        [vz_min,a,zp4_min]     = tracer_RK4(zp_min(:,N-1), ds, theta, bet, Q, ww, Qr);
        zp_min(:,N) = zp4_min + ds.*exp(i*(a));
        N = N+1;
	end
    
%Back transform the particle locations for the purpose of plotting them
        xp_t    = real(zp_maj);
        yp_t    = imag(zp_maj)/bet;
        xp      = xp_t * cos(-theta) + yp_t * sin(-theta);
        yp      = -xp_t * sin(-theta) + yp_t * cos(-theta);
        hold on;
        for jj  = 1:pnum
            plot(xp(jj,:),yp(jj,:),'k', 'LineWidth',2.5);
        end
        xp_t    = real(zp_min);
        yp_t    = imag(zp_min)/bet;
        xp      = xp_t * cos(-theta) + yp_t * sin(-theta);
        yp      = -xp_t * sin(-theta) + yp_t * cos(-theta);
        hold on;
        for jj  = 1:pnum
            plot(xp(jj,:),yp(jj,:),'k', 'LineWidth',2.5);
        end

%PLOT THE WELLS
        h = plot(wells(find(Qplot~=0),1),wells(find(Qplot~=0),2),'ko');
        set(h,'MarkerFaceColor','y');
        plot(wells(:,1),wells(:,2),'kx')
        h = plot(wells(find(Qplot==0),1),wells(find(Qplot==0),2),'ko');
        set(h,'MarkerFaceColor','r');     
%PLOT THE STAGNATION POINTS ON TOP SO THEY ARE CLEARLY VISIBLE
        h = plot(stag_x,stag_y,'kd');
        set(h,'MarkerFaceColor','g','MarkerSize',6);
% PLOT THE TERMINAL POINTS
%         x = real(z_stag_part_maj)*cos(-theta) + imag(z_stag_part_maj)/bet*sin(-theta);
%         y = -real(z_stag_part_maj)*sin(-theta) + imag(z_stag_part_maj)/bet*cos(-theta);
%         plot(x,y,'b.')
%         x = real(z_stag_part_min)*cos(-theta) + imag(z_stag_part_min)/bet*sin(-theta);
%         y = -real(z_stag_part_min)*sin(-theta) + imag(z_stag_part_min)/bet*cos(-theta);
%         plot(x,y,'r.')
        
axis([xmin xmax ymin ymax]);
axis equal
set(gca,'Color','None');
print(outfig,'-depsc',outfigname);

function [vz,a,zin4] = tracer_RK4(z_in, ds, theta, beta, Q, wells_Z_tr, Qrz)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  4th Order Runge-Kutta scheme used in ComCZAR_RK4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% streamline tracing by RK4 for tracer test

% predictor1
    W1 = tracer_discharge(z_in,wells_Z_tr, Q, Qrz);
    vz1 = conj(W1); 
    a   = angle(vz1);
    zin1    = z_in + ds/2.*exp(i*a);
% predictor2
    W2 = tracer_discharge(zin1,wells_Z_tr, Q, Qrz);
    vz2 = conj(W2); 
    a   = angle(vz2);
    zin2    = z_in + ds/2.*exp(i*a);
% precitor 3
    W3 = tracer_discharge(zin2,wells_Z_tr, Q, Qrz);
    vz3 = conj(W3);
    a   = angle(vz3);
    zin3    = z_in + ds.*exp(i*a);
% predictor 4
    W4 = tracer_discharge(zin3,wells_Z_tr, Q, Qrz);
    vz4 = conj(W4);
    % corrector
    vz = (vz1+2*vz2+2*vz3+vz4)/6;
    a = angle(vz); 
    zin4 = z_in;
    
function W= tracer_discharge(Z, wells_Z_tr, Q, Qrz)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to calculate complex discharge
% used by ComCZAR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

W = conj(Qrz);
for wellj = 1:length(wells_Z_tr)
    W = W - Q(wellj)/2/pi./(Z-wells_Z_tr(wellj)+eps);
end