function [] = profile51(fname,arcmap_directory,matlab_directory,exist_cd)

% profile51.m is the core stream profile analysis code. It works in 
% concert with sa_analysis51.m, sa_regress51.m, movavg51.m, usgscontour.m, 
% stepremover.m, answer_yn.m, slopes.m, closeto.m, closeto_allvalues.m,
%
% USAGE:
%       chandata = profile51('filename','arcmap_directory','matlab_directory','exist_cd');
%           exist_cd is a flag that indicates whether a new channel should
%           be extracted, or whether one intends to re-process
%           previously-extracted data
% INPUT:
%       'filename' is the part of the dem or acc file that does NOT contain -accm or -demm.
%           Because filename and directory_name are strings they MUST be written within single quotes.
%           It loads files exported as ASCII from ARC, headers trimmed, and saved as .mat files.
%           (these files are the masked dem and flowaccumulation matrices.  
%           Make these files manually, or use arcdemtxt2matlab.m).
%           filenameaccm.mat and filenamedemm.mat are loaded from the 'matlab_directory'.
%           The files must contain, respectively, variables named "accum" and "dem".
%                                           OR:  
%       'filename' is the prefix of an existing chandata file before _chandata.mat.
%       'arcmap_directory' and 'matlab_directory' are strings that give the complete directory (folder) path (see below).
%            Use '.' as shorthand to specify the local matlab directory.
%       NOTE: In addition to these input parameters, files
%       run_parameters.txt and location_ij.txt must exist in
%       'arcmap_directory'.  These files are created by the profile toolbar
%       in ArcGIS.
% OUTPUT:
%   chandata:  10 column matrix of DEM points along channel; columns are:
%        chandata = [dfd' pelev' drainarea' smooth_pelev' ptargi' ptargj' dfm' auto_ks_vals' x_coord' y_coord']
%        dfd = distance from divide
%        pelev = elevation
%        drainarea = drainage area
%        smooth_pelv = smoothed elevation
%        ptargi = x-coordinate (matrix value)
%        ptargj = y-coordinate (matrix value)
%        dfm = distance from mouth
%        auto_ks_vals = automatically extracted steepness indices along stream
%        x_coord = geographic coordinate (easting)
%        y_coord = geographic coordinate (northing)
%   Various matlab, postscript, and arcmap-loadable files are also interactively saved.
%
% 'arcmap_directory' is the directory where arcmap writes all stream-related text files and shapefiles.
% IT MUST MATCH THE DIRECTORY SPECIFIED IN ARCMAP ("set parameters to sent to matlab" button).
%
% Matlab writes files for input into arcmap into 'arcmap_directory' (e.g., name_stream.txt, section#.txt).
% 'matlab_directory' is the directory where new matlab files are saved.
 
% To turn on and off various interactive flags, search through profile51.m (and sa_analysis and sa_regress) to find "interactive parameter",
% and then comment out the question line and uncomment the appropriately-set
% value.  

% Commented OUT: EDITED to name stream text files as name_stream.txt, to
% work with 10-27-07 version of profiler.dll (not implemented at GSA Short Course)

% set global varible (tribsection).  This is used to number the
% files associated with the trib sections, which are read in Arcview
global TRIBSECTION

% set tribsection
TRIBSECTION = 1;

%flag to import geographic coordinates, rather than matlab indices:
use_geographic_coords=1;


%check number of input arguments; abort if not 6:
if exist_cd=='y',
    disp('Processing PREVIOUS chandata.')
    %       'exist_cd' is a flag to indicate whether a new chandata file is to be created (this is the default
    %           value: 'n', meaning there is no existing chandata file), or whether a pre-existing chandata file
    %           (with name fname_chandata) is to be processed ('y').
elseif exist_cd=='n',
    disp('Processing NEW chandata.')
else
    disp('Error--invalid entry for exist_cd, must be y or n.  Exiting')
    return
end
%expand directory strings if '.' was used:
if strcmp(arcmap_directory,'.'),
    arcmap_directory=pwd;
end
if strcmp(matlab_directory,'.'),
    matlab_directory=pwd;
end
%add a backslash to the directory strings if it isn't already there:
slash = '\';
if strcmp(arcmap_directory(length(arcmap_directory)),slash),
    arc_workdir = arcmap_directory;
else
    arc_workdir = strcat(arcmap_directory,slash);
end
if strcmp(matlab_directory(length(matlab_directory)),slash),
    mat_workdir = matlab_directory;
else
    mat_workdir = strcat(matlab_directory,slash);
end

%
% *************************************************************************
% Load Data Loop: Dem, Accum for exist_cd = 'n' (no existing chandata),
% Chandata for exist_cd = 'y'
% *************************************************************************
%

%loop to enter correct exist_cd param, if its wrong:
while exist_cd ~= 'n' & exist_cd ~= 'y', 
    exist_cd=input(sprintf('Use previous chandata file %s?  (y/n)',[fname,'_chandata;']),'s');
end

if exist_cd == 'n'
    % Load files saved to matlab binaries w/ matrices named dem, accum
    
    try
        eval(['load ',mat_workdir,fname,'demm;']);
    catch
        disp(sprintf('attempt to load dem file %s failed!',[mat_workdir,fname,'demm.mat;']));
        dem_fname = input(sprintf('Enter complete dem filename in %s:  ',mat_workdir),'s');
        eval(['load ',mat_workdir,dem_fname]);
    end
    
    try        
        eval(['load ',mat_workdir,fname,'accm']);
    catch
        disp(sprintf('attempt to load accum file %s failed!',[mat_workdir,fname,'accm.mat;']));
        acc_fname = input(sprintf('Enter complete accum filename in %s:  ',mat_workdir),'s');
        eval(['load ',mat_workdir,acc_fname]);
    end

    if use_geographic_coords,   %load *meta.mat file, that contains x,y--the geographic coordinates corresponding to the dem.
        try
            eval(['load ',mat_workdir,fname,'meta']);
            easting=x;  
            northing=y;
            clear x y 
        catch
            disp(sprintf('attempt to load metadata file %s failed!',[mat_workdir,fname,'meta.mat;']));
            meta_fname = input(sprintf('Enter complete metadata filename in %s:  ',mat_workdir),'s');
            eval(['load ',mat_workdir,meta_fname]);
            easting=x;  
            northing=y;
            clear x y 
        end
    end
        
    %error checking to make sure that variables named dem and accum now exist:  
    if isempty(who('accum')),
        disp(sprintf('Error--file %s apparently does not contain variable "accum".',[fname,'accm']));
%         disp('Manually create variable "accum", or load appropriate accumulation data, or type "dbquit" to exit profile3.m');
%         keyboard
        disp('Variable "accum" must exist; profile3.m exiting.')
        return
    end
    if isempty(who('dem')),
        disp(sprintf('Error--file %s apparently does not contain variable "dem".',[fname,'demm']));
%         disp('Manually create variable "dem", or load appropriate DEM data, or type "dbquit" to exit profile3.m')
%         keyboard
        disp('Variable "dem" must exist; profile3.m exiting.')
        return
    end    
    
    % load trib starting point, cellsize, and reference concavity from
    % arcview/arcgis
    
    eval(['load ', arc_workdir, 'run_parameters.txt']);
    cellsize = run_parameters(1);
    theta_ref = run_parameters(2);
    rmspike = run_parameters(3);
    no_step = run_parameters(4);
    smooth_prof = run_parameters(5);
    wind = run_parameters(6);
    cont_intv = run_parameters(7);
    ks_window = run_parameters(8);
    Pix_Ran_Dnst = run_parameters(9);
    MinAccum2start = run_parameters(10);
    gridsize = size(accum);
    
    movernset = -1*theta_ref;
    %gridsize used to test if targi,j stepping off grid
    
    % Set cellsize for analysis
    pix = round(10*cellsize)/10;
    diag = round(10*(sqrt((pix^2)+(pix^2))))/10;
    ar = round(pix^2);
    
elseif exist_cd == 'y'
    %
    % Else load existing chandata file; name = fname (for uniform input
    % convenience)
    %
    try
        eval(['load ',mat_workdir,fname,'_chandata;']);
    catch
        disp(sprintf('attempt to load chandata file %s failed!',[mat_workdir,fname,'_chandata.mat;']));
        chan_fname = input('Enter correct chandata filename:  ','s');
        try
            eval(['load ',mat_workdir,chan_fname]);
        catch
            disp(sprintf('Attempt to load chandata file %s failed again, exiting.',[mat_workdir,chan_fname]));
            return
        end
    end
    name = fname; 
    %
    %   read cellsize and reference concavity theta_ref from arcview directory
    %
    try
        eval(['load ', arc_workdir, 'run_parameters.txt'])
        cellsize = run_parameters(1);
        movernset = -1*run_parameters(2);
        rmspike = run_parameters(3);
        no_step = run_parameters(4);
        smooth_prof = run_parameters(5);
        wind = run_parameters(6);
        cont_intv = run_parameters(7);
        ks_window = run_parameters(8);
        Pix_Ran_Dnst = run_parameters(9);
        MinAccum2start = run_parameters(10);
    catch
        disp('Oops, can''t read file run_parameters.txt from arc_workdir.')
        cellsize=input('Enter cellsize in meters (usually 10 or 30):  ');
        movernset=input('Enter reference concavity (e.g. -0.4):   ');
    end

    pix = round(10*cellsize)/10;
    %
    
else
    text = 'No valid option specified. Input parameter exist_cd must be either n or y. Restart'
    return
end
%
% *************************************************************************
% End Load Data Loop
% *************************************************************************
%


%loop to load ALL points in location_ij.txt, FIRST TIME!  after first time, it can be re-read inside loop, see below.
bigloopnum=1;
if exist_cd == 'n',
        file = 'location_ij.txt';
        source = strcat(arc_workdir,file);
         %allow reading in and processing of multiple points
            % import format: [row column easting northing  name]
            % note: atargx = northing (rows), atargy = easting (columns)
        [atargiall,atargjall, atargyall,atargxall,nameall] = textread(source,'%f %f %f %f %s',-1);      %-1 flag will load all lines.  Assumes 5 columns
%         nameall = char(nameall);
        clear file source;
        numchanpts=size(atargiall,1); %number of rows in location_ij.txt, ie number of points to loop through
end


        
answer1 = 1;
while answer1,
    clear pdist paccum pelev drainarea targi targj ptargi ptargj dfm dfd p_x p_y
    %
    TRIBSECTION = 1;
    %
    % *************************************************************************
    % Create Chandata Loop for default, exist_cd = 'n' (no existing chandata)
    % *************************************************************************
    %
    if exist_cd == 'n'
        %
        disp('Computer is creating channel profile downstream of the point selected in ArcMap')

        if bigloopnum > numchanpts,   %case to reread location_ij.txt
            bigloopnum=1;
            file = 'location_ij.txt';
            source = strcat(arc_workdir,file);
            %allow reading in and processing of multiple points
            % import format: [row column easting northing  name]
            % note: atargx = northing (rows), atargy = easting (columns)
            [atargiall,atargjall,atargyall,atargxall,nameall] = textread(source,'%f %f %f %f %s',-1);      %-1 flag will load all lines.  Assumes 5 columns
%             nameall = char(nameall);
            clear file source;
            numchanpts=size(atargiall,1); %number of rows in location_ij.txt, ie number of points to loop through   
        end
        atargi=atargiall(bigloopnum);
        atargj=atargjall(bigloopnum);
        atargx=atargxall(bigloopnum);
        atargy=atargyall(bigloopnum);
        name=char(nameall(bigloopnum));
        disp(sprintf('Processing stream %s from ArcMap',name));

        if use_geographic_coords,       %convert location_ij just loaded from geographic coords to actually be matlab i,j
            %error check; make sure current point falls inside the bounds of the easting, northing vectors, ie in the bounds of the dem.
            if atargx<min(northing) | atargx>max(northing) | atargy<min(easting) | atargy>max(easting),
                disp('ERROR!  Starting point chosen, geographic coordinates, out of bounds of DEM.  Aborting.')
                warning('hi')
                keyboard
                return
            end
            northindex=closeto(atargx,northing);
            eastindex=closeto(atargy,easting);
            atargj=eastindex;
            atargi=northindex;
        end
        
        %
        targj=round(atargj);
        targi=round(atargi);
        
        % First, head DOWNSTREAM some # of steps (stepb)to make sure you are in the channel.
        for steps = 0:Pix_Ran_Dnst
            steps = steps+1;
            try
                patch = accum(targi-1:targi+1,targj-1:targj+1);
            catch
                warning('Error!  Most likely, channel starting point is not located on the DEM.  Aborting script.')
                return
            end
            [i_dnst,j_dnst] = find(patch==max(max(patch)));
            targi=targi+(i_dnst-2);
            targj=targj+(j_dnst-2);
        end

        % Now, search back UPSTREAM, following the path of largest upstream drainage areas.    
        upst_accum_present = accum(targi(1),targj(1));
        while upst_accum_present > MinAccum2start
            patch = accum(targi-1:targi+1,targj-1:targj+1); 
            list = sort(reshape(patch,9,1));
            [i_list,j_list] = find(list==upst_accum_present);
            output_val = list(i_list-1, j_list);
            [i_upst,j_upst] = find(patch==output_val(1));
            targi=targi+(i_upst-2);
            targj=targj+(j_upst-2);
            targi=targi(1);
            targj=targj(1);
            upst_accum_present = patch(i_upst(1),j_upst(1));
        end

        % Now work back DOWNSTREAM from the channel head to the outlet.
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        true=1;
		step=1;
        distance = [diag pix diag pix 0 pix diag pix diag];
        % Gather data from the first pixel of the channel
		pdist(step)=0;
		paccum(step)=accum(targi,targj);
		pelev(step)=dem(targi,targj);
%         plith(step)=rock(targi,targj);
		ptargi(step)=targi;
		ptargj(step)=targj;
    	
		% Begin the channel extraction while loop:
		while true == 1
			step=step+1;
            patch = accum(targi-1:targi+1,targj-1:targj+1);
            patch(2,2)=0;
            if step>2
                patch(2-(targi-ptargi(step-2)),2-(targj-ptargj(step-2)))=0;
            end
            [list,index] = sort(reshape(patch,9,1));
            dnst_accum = list(9);
            [i_dnst,j_dnst] = find(patch==dnst_accum);
            targi=targi+(i_dnst(1)-2);
            targj=targj+(j_dnst(1)-2);

            % now extract data for this new point
            pdist(step)=distance(index(9));
            paccum(step)=dnst_accum;
            pelev(step)=dem(targi,targj);
%             plith(step)=rock(targi,targj);
			ptargi(step)=targi;
			ptargj(step)=targj;
			p_x(step) = targi;
			p_y(step) = targj;
            
            % Note: will attempt to step out of bounds and will crash if last
            % point on channel is at edge of dem (i.e., if there is no NODATA
            % buffer)
            if paccum(step) < paccum(step-1)
                true = 0;
            elseif targi >= gridsize(1)
                true = 0;
            elseif targj >= gridsize(2)
                true = 0;
            elseif targi <=1
                true = 0;
            elseif targj <=1
                true = 0;
            end
    end
            	
        % Calculate cumulative distance from divide and convert to distance from the mouth.
        % Flip pelev, dfd  and paccum vectors.
        dfd = cumsum (pdist);
		pelev = fliplr(pelev);
%         plith = fliplr(plith);
		paccum = fliplr(paccum);
		dfd = fliplr(dfd);
		ptargi = fliplr(ptargi);
		ptargj = fliplr(ptargj);
		p_x = fliplr(p_x);
		p_y = fliplr(p_y);
		dfm = max(dfd)-dfd;
		% Convert drainage area from pix to m^2.  Create dummy data for smooth_pelev
		drainarea = (paccum.*ar);
		smooth_pelev = zeros(1,max(size(pelev)));  		
        auto_ks_vals = zeros(1,max(size(pelev)));   
		x_coord = zeros(1,max(size(pelev)));
		y_coord = zeros(1,max(size(pelev)));       
        

        % make the chandata matrix
        chandata = [dfd' pelev' drainarea' smooth_pelev' ptargi' ptargj' dfm' auto_ks_vals' x_coord' y_coord'];
    end
    %    
    % *************************************************************************
    % END Create Chandata Loop for default, exist_cd = 'n' (no existing chandata)
    % *************************************************************************
    %
    % *************************************************************************
    % Chandata now loaded for either option, begin analysis
    % *************************************************************************
    %

    
    % reset variables to ensure all are in the same row vector format.
    dfd = chandata(:,1)';
    pelev = chandata(:,2)';
    drainarea = chandata(:,3)';
    smooth_pelev = chandata(:,4)';
    ptargi = chandata(:,5)';
    ptargj = chandata(:,6)';
    p_x = ptargi;
    p_y = ptargj;
    dfm = chandata(:,7)';
    if exist_cd == 'n',
        if use_geographic_coords,    %convert j and i locations to x and y locations using the meta.mat file data
            chandata(:,9) = easting(ptargj,1);
            chandata(:,10) = northing(ptargi,1);
        end
    end
    % *************************************************************************
    % Loop to perform analysis starts here
    % *********************************************************************
    %

    if rmspike,
        temp_elev = fliplr(pelev);
        q = length(temp_elev)-1;
        for i = 1:q
            if temp_elev(i+1) > temp_elev(i)
                temp_elev(i+1) = temp_elev(i);
            end
        end
        new_pelev = fliplr(temp_elev);
    else
        new_pelev = pelev;
    end   

    if no_step,
        smooth_pelev = new_pelev;
        chandata(:,4) = smooth_pelev';
        wind = 0;
    else
        %only allow smoothing if step remover will not be used:  smooth raw
        %dem data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
        if smooth_prof,
            [smooth_pelev] = movavg51(new_pelev,pix,wind);
            % replace with option movavg51a to break into segments, avoid
            % smoothing across tributary junctions
            % [smooth_pelev] = movavg51a(new_pelev,drainarea,pix,wind);
            chandata(:,4) = smooth_pelev';
        else
            smooth_pelev = new_pelev;
            chandata(:,4) = smooth_pelev';
            wind = 0; %no smoothing done
        end
    end
    %
    % end spike removal / step removal / smoothing routine
    
    % write text file of necessary chandata out for Arcview.  This will be loaded into table for creating the stream data.
%     toGisData  = [dfd' pelev' drainarea' smooth_pelev' p_x' p_y']; %(the old format)
    toGisData  = [drainarea' chandata(:,8) p_x' p_y' chandata(:,9) chandata(:,10)];
    %
    [v,d]=version;
    if str2num(v(1))==6
        %dlmwrite([arc_workdir,name,'.txt'], toGisData,' ');
        dlmwrite([arc_workdir,name,'_stream.txt'], toGisData,' ');
    else
        %dlmwrite([arc_workdir,name,'.txt'],toGisData,'delimiter',' ','precision','%1.14g');
        dlmwrite([arc_workdir,name,'_stream.txt'],toGisData,'delimiter',' ','precision','%1.14g');
    end
%
    % Second: call function sa_analysis to run slope-area analysis and
    % knickpoint picking routine
    %

    [ida,islope,lbarea,lbslope,chandata,ans2] = sa_analysis51(chandata,movernset,name,arc_workdir,mat_workdir,rmspike,wind,no_step,ks_window,cont_intv);
    
    % *********************************************************************
    % BEGIN FINAL DATA/FIGURE-SAVING ROUTINE
    % *********************************************************************
 
%%%%%%%%%%%INTERACTIVE PARAMETER--Comment out interactive, uncomment set value, if desired    
    disp('Set the axis limits of the first two plots:')
    ans3 = answer_yn('Change the axis limits for fig 2, subplots 1 and 2?');
%    ans3 = 0;    %for no 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ans3,
        if ~isempty(who('min_elev')) & ~isempty(who('max_elev')) & ~isempty(who('max_length')),
            ans4=answer_yn('Use previous bounds?');
            if ans4,
            else
                min_elev=input('Enter minimum elevation (m):  ');
                max_elev=input('Enter maximum elevation (m):  ');
                max_length=input('Enter max channel length (km):  ');
            end
        else
            min_elev=input('Enter minimum elevation (m):   ');
            max_elev=input('Enter maximum elevation (m):   ');
            max_length=input('Enter max channel length (km):    ');
        end
        
        figure (2)
        subplot(3,1,1)
        try
            axis([0 max_length min_elev max_elev])
        catch
            disp('Error rescaling axis, probably from bad user input on min and max values.')
            disp('One more chance to enter reasonable values:')
            min_elev=input('Enter minimum elevation (m):   ');
            max_elev=input('Enter maximum elevation (m):   ');
            max_length=input('Enter max channel length (km):    ');
            try
                axis([0 max_length min_elev max_elev])
            catch
                disp('Error rescaling axis, probably from bad user input on mins and maxes.  Axes unchanged.')
            end
        end
        subplot(3,1,2)
        w2=axis;
        try
            axis([ 0 max_length w2(3) w2(4)])
        catch
            disp('Error rescaling axis, probably from bad user input on mins and maxes.  Axes unchanged.')
        end
       
        figure (3)
        subplot(3,1,1)
        w2=axis;
        try
            axis([0 max_length w2(3) w2(4)])
        catch
            disp('Error rescaling axis, probably from bad user input on mins and maxes.  Axes unchanged.')
        end
        subplot(3,1,2)
        w2=axis;
        try
            axis([0 max_length w2(3) w2(4)])
        catch
            disp('Error rescaling axis, probably from bad user input on mins and maxes.  Axes unchanged.')
        end
   
        subplot(3,1,3)
        w2=axis;
        try
            axis([0 max_length w2(3) w2(4)])
        catch    
            disp('Error rescaling axis, probably from bad user input on mins and maxes.  Axes unchanged.')
        end
    end
    
    %
    % Different save filename options depending on whether a pre-existing
    % chandata file is being reprocessed
    %
    if exist_cd == 'n'
%        t5 = answer_yn('Save this long profile (chandata.mat) and print figures?');
        t5 = ans2;
        if t5,
            sa_data = [ida' islope'];
            lb_sa_data = [lbarea' lbslope'];
            eval ([' save ',mat_workdir,name,'_chandata.mat chandata -MAT'])
            %        eval ([' save ',mat_workdir,name,'_chandata_ns.mat chandata_ns -MAT'])
            eval ([' save ',mat_workdir,name,'_sa_data.mat sa_data -MAT'])
            eval ([' save ',mat_workdir,name,'_lb_sa_data.mat lb_sa_data -MAT'])
        end
%        t3 = answer_yn('Print matlab figures (2 and 3) to file?');
        t3=ans2;
        if t3,
            figure(2)
            % arc_workdir,name,'.txt'  saveas(gcf,'pred_prey.fig')
            eval ([' print -depsc ',mat_workdir,name,'.eps'])
            figname = [arc_workdir name '.jpg'];
            saveas (gcf,figname)
            figure(3)
            eval ([' print -depsc ',mat_workdir,name,'_f3.eps'])
            figname = [arc_workdir name '_f3.jpg'];
            saveas (gcf,figname)
        end
        
        disp('You MUST return to ArcMap and import data using the "add stream shapefile"')
        disp('and "create knickpoint shapefile" buttons.  Otherwise NO SHAPEFILES will be saved.')
        bigloopnum=bigloopnum+1;
        if bigloopnum <= numchanpts,    %case where there are more rows in location_ij.txt to read as channel heads
            disp('After importing data to ArcMap, the next point channel starting point previously chosen will be used (from location_ij.txt).')
%             answer1 = answer_yn(sprintf('Continue with the next starting point, %s, previously chosen from ArcMap? (y to continue, n to quit)',nameall(bigloopnum,:)));
            answer1 = answer_yn(sprintf('Continue with the next starting point, %s, previously chosen from ArcMap? (y to continue, n to quit)',char(nameall(bigloopnum))));
        else
            disp('No more points chosen for processing!')
            disp('After importing data, you MUST select another tributary starting point ("send a point to matlab" button) in ArcMap.')
            disp('Otherwise, previous point will be used.')
            answer1 = answer_yn('Are you ready to grab another starting point from ArcMap? (y to continue, n to quit)');
        end
        
    else if exist_cd == 'y'
            
            suffix = '';
            t5 = answer_yn('Save this processed long profile matlab data?');
            if t5,
                suffix = input(sprintf('Enter optional suffix (*) for naming %s_*_sa_data, %s_*_lb_sa_data and %s_*_chandata :  ',name,name,name),'s'); 
                sa_data = [ida' islope'];
                lb_sa_data = [lbarea' lbslope'];
                eval ([' save ',mat_workdir,name,'_',suffix,'_chandata.mat chandata -MAT'])
                %%right now saving doesn't save all of the no-step or interpolated chandata points; only saves the slope and area at these points.
                %        eval ([' save ',mat_workdir,name,'_chandata_ns.mat chandata_ns -MAT'])
                %sa_data is 2 columns, 1st col drainage area, 2nd column slope, only at the chosen (step removed or interpolated) subset of points.
                eval ([' save ',mat_workdir,name,'_',suffix,'_sa_data.mat sa_data -MAT'])
                %lb is for log-binned, 2 col drainage area and slope (as shown in slope-area plot, fig2 plot3).
                eval ([' save ',mat_workdir,name,'_',suffix,'_lb_sa_data.mat lb_sa_data -MAT'])
            end
            t3 = answer_yn('Print matlab figures (2 and 3) to file?');
            if t3,
                figure(2)
                eval ([' print -depsc ',mat_workdir,name,'_',suffix,'_fig2.ps'])
                figure(3)
                eval ([' print -depsc ',mat_workdir,name,'_',suffix,'_fig3.ps'])
            end
            %
            % no repeat cycle if reprocessing specified chandata file
            %
            answer1 = 0;
            
        end
        %
        % End Save Options Loop
        %
    end
    %
    % End While Loop for repeat analysis if answer1 = 1, END profile3.m
    %
end