clear

 vhe_model='calc01';    % Tsim input data model to be used in this wear prediction loop
 vhe_model='calc01';    % R600 stiff bogie
%vhe_model='calc02';    % Tangent track; No track irr.; With traction
 iloop=  0;             % Start wheel profile #
 nloop=  3;             % End wheel profile #


 string= sprintf('mkdir -p qout');
 [istat, output]= system(string); if (istat ~= 0); error(output); end

 string= sprintf('mkdir -p kpf/qout');
 [istat, output]= system(string); if (istat ~= 0); error(output); end

 ident= sprintf('%s_%3.3d', vhe_model, iloop);

 if (iloop==0)
  string= ['cp -p kpf/w_prof/outer_wheel_000.wheel kpf/w_prof/outer_wheel_', ident, '.wheel'];
  [istat, output]= system(string); if (istat ~= 0); error(output); end
  string= ['cp -p kpf/w_prof/inner_wheel_000.wheel kpf/w_prof/inner_wheel_', ident, '.wheel'];
  [istat, output]= system(string); if (istat ~= 0); error(output); end
 end

 outer_wheel_prof= sprintf('outer_wheel_%s', ident);
 inner_wheel_prof= sprintf('inner_wheel_%s', ident);
 string= ['[outer_wheel_prof_X, outer_wheel_prof_Y]= read_2col_file ("kpf/w_prof/',outer_wheel_prof,'.wheel");'];   % Read the initial profiles
 try
     eval(string)
 catch
     error("Exit status from command %s\n       not equal zero\n%s\n", string, lasterr ());
 end
 string= ['[inner_wheel_prof_X, inner_wheel_prof_Y]= read_2col_file ("kpf/w_prof/',inner_wheel_prof,'.wheel");'];
 try
     eval(string)
 catch
     error("Exit status from command %s\n       not equal zero\n%s\n", string, lasterr ());
 end

%
%  Wear step loop
%  ------------------------
while( iloop<nloop ) % {

%
disp('Run kpfr     {{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{')
 ident= sprintf('%s_%3.3d', vhe_model, iloop);
 outer_wheel_prof= sprintf('outer_wheel_%s', ident);
 inner_wheel_prof= sprintf('inner_wheel_%s', ident);

 cd kpf
  tic;
  string= ['sed -e "s!WHEEL!w_prof/',outer_wheel_prof,'.wheel!" -e "s!RAIL!$genkpf/../r_prof/uic60/uic60i40.rail!"  kpff/Master_same_profiles.kpff > kpff/',outer_wheel_prof,'.kpff'];
  [istat, output]= system(string); if (istat ~= 0); error(output); end
  string= ['kpf kpff/',outer_wheel_prof,'.kpff > qout/',outer_wheel_prof,'.out'];
  [istat, output]= system(string); if (istat ~= 0); error(output); end
  string= ['sed -e "s!WHEEL!w_prof/',inner_wheel_prof,'.wheel!" -e "s!RAIL!$genkpf/../r_prof/uic60/uic60i40.rail!"  kpff/Master_same_profiles.kpff > kpff/',inner_wheel_prof,'.kpff'];
  [istat, output]= system(string); if (istat ~= 0); error(output); end
  string= ['kpf kpff/',inner_wheel_prof,'.kpff > qout/',inner_wheel_prof,'.out'];
  [istat, output]= system(string); if (istat ~= 0); error(output); end
 cd ..
 t_kpfr= toc;
disp(['t_kpfr= ',num2str(round(t_kpfr)),' s'])
disp('End kpfr     }}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}')



%
disp('Run tsim     {{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{');
  tic;
  string= sprintf('sed "/Wheel-rail geometry functions/ a\\\n');
  strin2= sprintf('insert file kpf/kpfr/%s.kpfr \\\\n', outer_wheel_prof);
  string= strcat(string,strin2);
  strin2= sprintf('insert file kpf/kpfr/%s.kpfr \\\\n', inner_wheel_prof);
  string= strcat(string,strin2);
  strin2= sprintf('in_substruct %s [ 111 ] \\\\n', outer_wheel_prof);
  string= strcat(string,strin2);
  strin2= sprintf('in_substruct %s [ 112 ] \\\\n', inner_wheel_prof);
  string= strcat(string,strin2);
  strin2= sprintf('in_substruct %s [ 121 ] \\\\n', inner_wheel_prof);
  string= strcat(string,strin2);
  strin2= sprintf('in_substruct %s [ 122 ]"',      outer_wheel_prof);
  string= strcat(string,strin2);
  strin2= sprintf(' runf/%s.tsimf > runf/%s.tsimf',  vhe_model, ident);
  string= strcat(string,strin2);

  [istat, output]= system(string); if (istat ~= 0); error(output); end

  string= sprintf('tsim -no_addarg -overwrite runf/%s.tsimf > qout/tsim_%s.out', ...
                                                  ident,      ident);
  [istat, output]= system(string); if (istat ~= 0); error(output); end
  t_tsim= toc;
disp(['t_tsim= ',num2str(round(t_tsim)),' s'])
disp('End tsim     }}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}')



%
disp('Run mplot    {{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{')
  tic;
  [istat, output]= system('del -q mplotr/wear_111l.print'); if (istat ~= 0); error(output); end
  [istat, output]= system('del -q mplotr/wear_111r.print'); if (istat ~= 0); error(output); end
  [istat, output]= system('del -q mplotr/wear_112l.print'); if (istat ~= 0); error(output); end
  [istat, output]= system('del -q mplotr/wear_112r.print'); if (istat ~= 0); error(output); end
  [istat, output]= system('del -q mplotr/wear_121l.print'); if (istat ~= 0); error(output); end
  [istat, output]= system('del -q mplotr/wear_121r.print'); if (istat ~= 0); error(output); end
  [istat, output]= system('del -q mplotr/wear_122l.print'); if (istat ~= 0); error(output); end
  [istat, output]= system('del -q mplotr/wear_122r.print'); if (istat ~= 0); error(output); end

  string= sprintf('mplot -no_addarg -no_interactive -no_use_MPfile mplotf/wear.mplotf %s > qout/mplot_%s.out', ident, ident);
  [istat, output]= system(string); if (istat ~= 0); error(output); end
  t_mplot=toc;
disp(['t_mplot= ',num2str(round(t_mplot)),' s'])
disp('End mplot    }}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}')


%
disp('Update wheel profile  {{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{')

  load -ascii mplotr/wear_111l.print
  load -ascii mplotr/wear_111r.print
  load -ascii mplotr/wear_112l.print
  load -ascii mplotr/wear_112r.print
  load -ascii mplotr/wear_121l.print
  load -ascii mplotr/wear_121r.print
  load -ascii mplotr/wear_122l.print
  load -ascii mplotr/wear_122r.print

% wear_111l = importdata('mplotr/wear_111l.print', ' ', 12);
% wear_111r = importdata('mplotr/wear_111r.print', ' ', 12);
% wear_112l = importdata('mplotr/wear_112l.print', ' ', 12);
% wear_112r = importdata('mplotr/wear_112r.print', ' ', 12);
% wear_121l = importdata('mplotr/wear_121l.print', ' ', 12);
% wear_121r = importdata('mplotr/wear_121r.print', ' ', 12);
% wear_122l = importdata('mplotr/wear_122l.print', ' ', 12);
% wear_122r = importdata('mplotr/wear_122r.print', ' ', 12);

% wear_111l = wear_111l.data;
% wear_111r = wear_111r.data;
% wear_112l = wear_112l.data;
% wear_112r = wear_112r.data;
% wear_121l = wear_121l.data;
% wear_121r = wear_121r.data;
% wear_122l = wear_122l.data;
% wear_122r = wear_122r.data;

  wear_outer= (wear_111l(:,2) + wear_111r(:,2) + wear_122l(:,2) + wear_122r(:,2)) /4.*1000.*100.;  % Convert wear to mm
  wear_inner= (wear_112l(:,2) + wear_112r(:,2) + wear_121l(:,2) + wear_121r(:,2)) /4.*1000.*100.;  % Multiply with 100 in order to
                                                                                                   % fake a longer simulation
  outer_wheel_prof_Y= outer_wheel_prof_Y - wear_outer;  %   update profile
  inner_wheel_prof_Y= inner_wheel_prof_Y - wear_inner;

  ident= sprintf('%s_%3.3d', vhe_model, iloop+1);
  write_2col_file(['kpf/w_prof/outer_wheel_', ident, '.wheel'], outer_wheel_prof_X, outer_wheel_prof_Y)
  write_2col_file(['kpf/w_prof/inner_wheel_', ident, '.wheel'], inner_wheel_prof_X, inner_wheel_prof_Y)

disp('End    wheel profile  }}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}')

  iloop= iloop + 1;
 end % } ===================================================================================


 figure
 hold on
 plot(outer_wheel_prof_X,outer_wheel_prof_Y)
 hold off
