function  combine_on_y (file_1, file_2, file_3, file_4)
% Copyright 2000, J.E. Akin.  All rights reserved.
%  ------------------------------------------------------
%   combine Cosine hill results for Paper
%  ------------------------------------------------------

% c_x    = x coordinates of nod_per_el line polygon
% c_y    = y coordinates of nod_per_el line polygon
% msh_typ_nodes   = connectivity list for elements, nt x nod_per_el
% nod_per_el = Nodes per element
% np     = Number of Points
% nt     = Number of elements
% pre_e  = Element items before connectivity list
  pre_e  = 0 ;
% pre_p  = Nodal items before coordinates
  pre_p = 1;
% msh_bc_xyz = Nodal coordinates (with preceeding data)
% t_x    = x coordinates of nod_per_el corners
% t_y    = y coordinates of nod_per_el corners
  BOG = -8 ; % 1st of graph tag to sepreate file saves
  EOG = -9 ; % end of graph tag to sepreate file saves

format short g
clear fid_1 fid_2 fid_3 fid_4

fid_X = fopen('Q_x_coord','r') ; % fudge incorrect x values
Q_1 = fscanf (fid_X, '%g %g %g \n', [3 inf]) ;
Q_1 = Q_1' ; % transpose back to rows
 
  if ( nargin == 1 )
    fid_1 = fopen(file_1, 'r') ;
    fid_2 = 0 ;
    fid_3 = 0 ;
    fid_4 = 0 ;
    file_2 = ['null'] ; 
    file_3 = ['null'] ;
    file_4 = ['null'] ;
  elseif ( nargin == 2 )
    fid_1 = fopen(file_1, 'r') ;
    fid_2 = fopen(file_2, 'r') ;
    fid_3 = 0 ;
    fid_4 = 0 ;
    file_3 = ['null'] ;
    file_4 = ['null'] ;
  elseif ( nargin == 3 )
    fid_1 = fopen(file_1, 'r') ;
    fid_2 = fopen(file_2, 'r') ;
    fid_3 = fopen(file_3, 'r') ;
    fid_4 = 0 ;
    file_4 = ['null'] ;
  else
    fid_1 = fopen(file_1, 'r') ;
    fid_2 = fopen(file_2, 'r') ;
    fid_3 = fopen(file_3, 'r') ;
    fid_4 = fopen(file_4, 'r') ;
  end % if no arguments
% display(file_1) ; display( file_2) ; 
% display( file_3) ; display( file_4)

% fid_2, fid_3, fid_4
%    load the file array   
 T_1 = fscanf (fid_1, '%g %g %g \n', [3 inf]) ;
 T_1 = T_1' ; % transpose back to rows

%    Set control data: number of points
  nr = size (T_1, 1) ;    % number of nodal points  
  %b fprintf ('Read %g triplets \n', nr)
%  display (T_1(1, 1:3))  % start_flag, component, x
%  display (T_1(2, 1:3))  % yr mon day saved
%  display (T_1(3, 1:3))  % hr min sec saved
%b display (T_1(4, 1:3))  % nt np per mesh   
%  display (T_1(nr, 1:3))  % end_flag min_v max_v
  nt  = 576 ; np = 625 ;   % hard code size of mesh
  V_N = T_1 (nr, 2) ; V_X = T_1 (nr, 3) ;
  nd = nr - 4 ; % data pairs  %b --------- 5 ---- later
  
%    Strip out the y-coordinates and values
  X_1 = T_1 (1, 3)  ; % x where saved
  %b y_1 (1:nd) = T_1 ((1:nd)+3, 2) ; % coord  %b --------- 5 ---- later
  y_1 (1:nd) = Q_1 ((1:nd)+3, 2) ; % coord  %b --------- 5 ---- later
  v_1 (1:nd) = T_1 ((1:nd)+3, 3) ; % coord  %b --------- 5 ---- later
  y_min = min (y_1) ;
  y_max = max (y_1) ;

  if ( fid_2 > 0 )
    T_2 = fscanf (fid_2, '%g %g %g \n', [3 inf]) ;
    T_2 = T_2' ; % transpose back to rows
    v_2 (1:nd) = T_2 ((1:nd)+3, 3) ;
    if ( T_2 (nr, 2) < V_N )
      V_N = T_2 (nr, 2) ;
    end % if
    if ( T_2 (nr, 3) > V_X )
      V_X = T_2 (nr, 3) ;
    end % if
  end % if

  if ( fid_3 > 0 )
    T_3 = fscanf (fid_3, '%g %g %g \n', [3 inf]) ;
    T_3 = T_3' ; % transpose back to rows
    v_3 (1:nd) = T_3 ((1:nd)+3, 3) ;
    if ( T_3 (nr, 2) < V_N )
      V_N = T_3 (nr, 2) ;
    end % if
    if ( T_3 (nr, 3) > V_X )
      V_X = T_3 (nr, 3) ;
    end % if
  end % if

  if ( fid_4 > 0 )
    T_4 = fscanf (fid_4, '%g %g %g \n', [3 inf]) ;
    T_4 = T_4' ; % transpose back to rows
    v_4 (1:nd) = T_4 ((1:nd)+3, 3) ;
    if ( T_4 (nr, 2) < V_N )
      V_N = T_4 (nr, 2) ;
    end % if
    if ( T_4 (nr, 3) > V_X )
      V_X = T_4 (nr, 3) ;
    end % if
  end % if

%    Get the exact values
  dy  = (y_max - y_min)/ 200  ;
  y_e = [y_min:dy:y_max] ;
  for j = 1:size (y_e, 2) ;
    r = sqrt (X_1 ^2 + y_e (j) ^2) ;
    if ( r < 0.5 )
      v_e (j) = sin ( pi*(1 - 2*r) ) ;
    else
      v_e (j) = 0 ;
    end % if
  end % for
  e_min = min (v_e) ; e_max = max (v_e) ;
  if ( e_min < V_N )
    V_N = e_min ;
  end % if
  if ( e_max > V_X )
    V_X = e_max ;
  end % if

%  set up plots
  clf

  plot(y_e, v_e, 'k-')
  hold on                           % hold image for plots
  grid
  plot(y_1, v_1, 'ko') % . o x + * s d v ^ < > p h
  if ( fid_2 > 0 )
    plot(y_1, v_2, 'ks') % . o x + * s d v ^ < > p h
  end % if
  if ( fid_3 > 0 )
    plot(y_1, v_3, 'k+') % . o x + * s d v ^ < > p h
  end % if
  if ( fid_4 > 0 )
    plot(y_1, v_4, 'k<') % . o x + * s d v ^ < > p h
  end % if

  axis ([y_min y_max V_N V_X])
  %b axis ([y_min y_max 0.0 1.1])

  %b xlabel (['X for points at Y = ', num2str(X_1)])      
  %b ylabel ('Solution Values')
  %b title(['Expected and Nodal FE Solutions', ...     
         %b '  (', int2str(nt),' Elements, ', ...
         %b int2str(np),' Nodes)'])

  if ( nargin == 1)
    legend ('Exact', 'Element S1 Q16  ' )
  elseif ( nargin == 2)
  % legend ('Exact', 'Element S1  Q16' , 'Element UGN Q16  ' )
  % legend ('Exact', 'Element SA2 Q4  ', 'Element SA2 Q9  ' )
  % legend ('Exact', 'Nodal SA2 Q4  ', 'Nodal SA2 Q9  ' )
  % legend ('Exact', 'Nodal S1 T3  ', 'Nodal SA2 T3  ' )
  % legend ('Exact', 'Nodal S1 T6  ', 'Nodal SA2 T6  ' )
    %legend ('Exact', 'Nodal S1 T10 ', 'Nodal SA2 T10 ' )
    %legend ('Exact', 'S1 QP T3', 'S1 QP Q4')
    %legend ('Exact', 'S1 QP T6', 'S1 QP Q9')
    legend ('Exact', 'S1 QP T10', 'S1 QP Q16')
  elseif ( nargin == 3)
  % legend ('Exact', 'Element S1  Q16' , 'Element UGN Q16' , ...
  %         'Nodal S1  Q16  ' )
  % legend ('Exact', 'Element S1 Q4  ' , 'Element S1 Q9  ' , ...
  %         'Element S1 Q16  ' )
  % legend ('Exact', 'Nodal S1 Q4  ' , 'Nodal S1 Q9  ' , ...
  %         'Nodal S1 Q16  ' )
  % legend ('Exact', 'Element UGN Q4  ', 'Element UGN Q9  ', ...
  %         'Element UGN Q16  ' )
  % legend ('Exact', 'Element S1  T3  ', 'Element UGN T3  ', ...
  %         'Element SA2 T3  ' )
  % legend ('Exact', 'Element S1  T6  ', 'Element UGN T6  ', ...
  %         'Element SA2 T6  ' )
  % legend ('Exact', 'S1 QP T3', 'S1 QP T6', 'S1 QP T10') 
  % legend ('Exact', 'S1 QP Q4', 'S1 QP Q9', 'S1 QP Q16' )
  % legend ('Exact', 'Q4 SUGN1', 'Q4 S1 E', 'Q4 S1 QP')    
    legend ('Exact', 'T3 SUGN1', 'T3 S1 E', 'T3 S1 QP')    
  elseif ( nargin == 4)
  % legend ('Exact', 'Element S1  Q16' , 'Element UGN Q16' , ...
  %         'Nodal S1  Q16' , 'Nodal SA2 Q16  ' )
  % legend ('Exact', 'Element S1  T10' , 'Element UGN T10' , ...
  %         'Nodal S1  T10' , 'Nodal SA2 T10  ' )
  % legend ('Exact', 'Element S1  T6' , 'Element UGN T6' , ...
  %         'Nodal S1  T6' , 'Nodal SA2 T6  ' )
    legend ('Exact', 'S1 QP Q4', 'S1 QP Q8', 'S1 QP Q9',   ...
            'S1 QP Q16' )
  end % if

%       -depsc -tiff            % for an eps version
  %b print ('-dpsc', ['combined_', int2str(i_p), '_on_y'])
  hold off
  %b v_text = ['Created combined_', int2str(i_p), '_on_y.ps'] ;
  %b fprintf (1,'%s', v_text) ; fprintf (1, ' \n' )
% end of combine_on_y