clear

% octave -q --no-gui
% wear_loop

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

 string= sprintf('mkdir -p kpf/qout');
 [istat, output]= system(string); if (istat ~= 0); error(output); exit; 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); exit; 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); exit; 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); exit; end
%
  string= ['kpf kpff/',outer_wheel_prof,'.kpff > qout/',outer_wheel_prof,'.out'];               %% Run program KPF outer_wheel_prof
  [istat, output]= system(string);
  if (istat ~= 0)
   printf(' ***ERROR*** In wear_loop.m\n');
   printf('             Exit status from program KPF not equal 0 (zero)\n');
   printf('             Please open and read the error message in file:\n');
   printf('             kpf/qout/%s.out\n', outer_wheel_prof);
   exit
  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); exit; end
%
  string= ['kpf kpff/',inner_wheel_prof,'.kpff > qout/',inner_wheel_prof,'.out'];               %% Run program KPF inner_wheel_prof
  [istat, output]= system(string);
  if (istat ~= 0)
   printf(' ***ERROR*** In wear_loop.m\n');
   printf('             Exit status from program KPF not equal 0 (zero)\n');
   printf('             Please open and read the error message in file:\n');
   printf('             kpf/qout/%s.out\n', inner_wheel_prof);
   exit
  end
 cd ..
 t_kpfr= toc;
disp(['t_kpfr= ',num2str(round(t_kpfr)),' s'])
disp('End kpfr     }}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}')


% Adjust the wear depending on contact angle
disp('Beg gamma_cos }}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}')

  string= ['(cd kpf/kpfr; kpf_conv_kpfr2matlab ',outer_wheel_prof,'.kpfr ',outer_wheel_prof,'.m)']
  [istat, output]= system(string);
  if (istat ~= 0)
   printf(' ***ERROR*** In wear_loop.m\n');
   printf('             Exit status from program "kpf_conv_kpfr2matlab" not equal 0 (zero)\n');
   printf('             output= %s\n', output);
   exit
  end
  eta= [-.040:.001:.040];
  cd kpf/kpfr
  eval (outer_wheel_prof)
  cd ../..

  cp1_zfn= interp1 (cp1_zfn_(:,1),cp1_zfn_(:,2), eta, "extrap");
  cp2_zfn= interp1 (cp2_zfn_(:,1),cp2_zfn_(:,2), eta, "extrap");
  cp3_zfn= interp1 (cp3_zfn_(:,1),cp3_zfn_(:,2), eta, "extrap");
  cp1_gamfn = interp1 (cp1_gamfn_(:,1),cp1_gamfn_(:,2), eta, "extrap");
  cp2_gamfn = interp1 (cp2_gamfn_(:,1),cp2_gamfn_(:,2), eta, "extrap");
  cp3_gamfn = interp1 (cp3_gamfn_(:,1),cp3_gamfn_(:,2), eta, "extrap");
  cp1_poswfn= interp1 (cp1_poswfn_(:,1),cp1_poswfn_(:,2), eta, "extrap");
  cp2_poswfn= interp1 (cp2_poswfn_(:,1),cp2_poswfn_(:,2), eta, "extrap");
  cp3_poswfn= interp1 (cp3_poswfn_(:,1),cp3_poswfn_(:,2), eta, "extrap");

% plot(eta,cp1_zfn)
% axis([ -.2 .2 -.005 .005 ])

  for ii= 1:length(eta)
   cpa_gamfn(ii) = cp1_gamfn(ii);
   cpa_poswfn(ii)= cp1_poswfn(ii);
   if (cp2_zfn(ii) < cp1_zfn(ii))
    cpa_gamfn(ii) = cp2_gamfn(ii);
    cpa_poswfn(ii)= cp2_poswfn(ii);
   endif
   if (cp3_zfn(ii) < cp2_zfn(ii))
    cpa_gamfn(ii) = cp3_gamfn(ii);
    cpa_poswfn(ii)= cp3_poswfn(ii);
   endif
  endfor
  
% plot(eta,cpa_gamfn,';cpa;', eta,cp1_gamfn,';cp1;', eta,cp2_gamfn,';cp2;', eta,cp3_gamfn,';cp3;')
% plot(eta,cp1_gamfn,';cp1;', eta,cp2_gamfn,';cp2;', eta,cp3_gamfn,';cp3;', eta,cpa_gamfn,';cpa;')
% plot(eta,cp1_poswfn,';cp1;', eta,cp2_poswfn,';cp2;', eta,cp3_poswfn,';cp3;', eta,cpa_poswfn,';cpa;')
% plot(eta,cpa_poswfn)

  for ii= 2:length(eta)                         ;% Avoid backsteps
   if (cpa_poswfn(ii) <= cpa_poswfn(ii-1))
    cpa_poswfn(ii)= cpa_poswfn(ii-1)+eps;
   endif
  endfor

  posw_111l_1D= [ -4.475000E-02
-4.425000E-02
-4.375000E-02
-4.325000E-02
-4.275000E-02
-4.225000E-02
-4.175000E-02
-4.125000E-02
-4.075000E-02
-4.025000E-02
-3.975000E-02
-3.925000E-02
-3.875000E-02
-3.825000E-02
-3.775000E-02
-3.725000E-02
-3.675000E-02
-3.625000E-02
-3.575000E-02
-3.525000E-02
-3.475000E-02
-3.425000E-02
-3.375000E-02
-3.325000E-02
-3.275000E-02
-3.225000E-02
-3.175000E-02
-3.125000E-02
-3.075000E-02
-3.025000E-02
-2.975000E-02
-2.925000E-02
-2.875000E-02
-2.825000E-02
-2.775000E-02
-2.725000E-02
-2.675000E-02
-2.625000E-02
-2.575000E-02
-2.525000E-02
-2.475000E-02
-2.425000E-02
-2.375000E-02
-2.325000E-02
-2.275000E-02
-2.225000E-02
-2.175000E-02
-2.125000E-02
-2.075000E-02
-2.025000E-02
-1.975000E-02
-1.925000E-02
-1.875000E-02
-1.825000E-02
-1.775000E-02
-1.725000E-02
-1.675000E-02
-1.625000E-02
-1.575000E-02
-1.525000E-02
-1.475000E-02
-1.425000E-02
-1.375000E-02
-1.325000E-02
-1.275000E-02
-1.225000E-02
-1.175000E-02
-1.125000E-02
-1.075000E-02
-1.025000E-02
-9.750001E-03
-9.250000E-03
-8.749999E-03
-8.250002E-03
-7.750001E-03
-7.250000E-03
-6.749999E-03
-6.250001E-03
-5.750000E-03
-5.249999E-03
-4.749998E-03
-4.250001E-03
-3.750000E-03
-3.249999E-03
-2.749998E-03
-2.250001E-03
-1.750000E-03
-1.249999E-03
-7.500015E-04
-2.500005E-04
 2.500005E-04
 7.500015E-04
 1.249999E-03
 1.750000E-03
 2.250001E-03
 2.750002E-03
 3.249999E-03
 3.750000E-03
 4.250001E-03
 4.750002E-03
 5.249999E-03
 5.750000E-03
 6.250001E-03
 6.750003E-03
 7.250000E-03
 7.750001E-03
 8.250002E-03
 8.749999E-03
 9.250000E-03
 9.750001E-03
 1.025000E-02
 1.075000E-02
 1.125000E-02
 1.175000E-02
 1.225000E-02
 1.275000E-02
 1.325000E-02
 1.375000E-02
 1.425000E-02
 1.475000E-02
 1.525000E-02
 1.575000E-02
 1.625000E-02
 1.675000E-02
 1.725000E-02
 1.775000E-02
 1.825000E-02
 1.875000E-02
 1.925000E-02
 1.975000E-02
 2.025000E-02
 2.075000E-02
 2.125000E-02
 2.175000E-02
 2.225000E-02
 2.275001E-02
 2.325000E-02
 2.375000E-02
 2.425000E-02
 2.475000E-02
 2.525000E-02
 2.575000E-02
 2.625000E-02
 2.675000E-02
 2.725000E-02
 2.775000E-02
 2.825000E-02
 2.875000E-02
 2.925000E-02
 2.975000E-02
 3.025001E-02
 3.075000E-02
 3.125000E-02
 3.175000E-02
 3.225000E-02
 3.275000E-02
 3.325000E-02
 3.375001E-02
 3.425000E-02
 3.475000E-02
 3.525000E-02
 3.575000E-02
 3.625000E-02
 3.675000E-02
 3.725000E-02
 3.775001E-02
 3.825000E-02
 3.875000E-02
 3.925000E-02
 3.975000E-02
 4.025000E-02
 4.075000E-02
 4.125001E-02
 4.175000E-02
 4.225000E-02
 4.275000E-02
 4.325000E-02
 4.375000E-02
 4.425000E-02
 4.475001E-02];
  posw_eta= interp1 (cpa_poswfn,eta, posw_111l_1D, "extrap");
% for ii= 1:length(posw_111l_1D)
%  if (posw_eta(ii) < -.040) posw_eta(ii)= -.040; endif
%  if (posw_eta(ii) >  .040) posw_eta(ii)=  .040; endif
% endfor

  gamma= interp1 (eta,cpa_gamfn, posw_eta, "extrap");   ;% Pick up gamma
  for ii= 1:length(eta)                                 ;% Avoid too big angles in gamma
   if (gamma(ii) > 1.5)
    cpa_poswfn(ii)= 1.5;
   endif
  endfor
  gamma_cos= cos(gamma);
% plot(posw_eta,gamma)
% plot(posw_eta,gamma_cos)
% axis([ -40 40 -1 1 ])
disp('End gamma_cos }}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}')


%
disp('Run tsim     {{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{');
  tic;

# awk '/jack/ {c+=1} {if(c==2){sub("jack","jill",$0)};print}' file
# awk '/Wheel-rail geometry functions/ {c+=1} { if (c==2) {
# printf("insert file kpf/kpfr/%s.kpfr \\\\n", outer_wheel_prof);
# }}' file

  string= sprintf("awk '/Wheel-rail geometry functions/ {c+=1} {print; if (c==2) {++c;");
  strin2= sprintf("printf(\"insert file kpf/kpfr/%s.kpfr\\n\");", outer_wheel_prof);  string= strcat(string,strin2);
  strin2= sprintf("printf(\"insert file kpf/kpfr/%s.kpfr\\n\");", inner_wheel_prof);  string= strcat(string,strin2);
  strin2= sprintf("printf(\"in_substruct %s [ 111 ] \\n\");",     outer_wheel_prof);  string= strcat(string,strin2);
  strin2= sprintf("printf(\"in_substruct %s [ 112 ] \\n\");",     inner_wheel_prof);  string= strcat(string,strin2);
  strin2= sprintf("printf(\"in_substruct %s [ 121 ] \\n\");",     inner_wheel_prof);  string= strcat(string,strin2);
  strin2= sprintf("printf(\"in_substruct %s [ 122 ] \\n\")}}'",   outer_wheel_prof);  string= strcat(string,strin2);
  strin2= sprintf(' runf/%s.tsimf > runf/%s.tsimf',     vhe_model, ident);  string= strcat(string,strin2)

# strin2= sprintf("printf(\"runf/%s.tsimf > runf/%s.tsimf",     vhe_model, ident);  string= strcat(string,strin2);
# strin2= sprintf("printf(\"insert file kpf/kpfr/%s.kpfr \\\\n\", inner_wheel_prof)");  string= strcat(string,strin2);
# strin2= sprintf(' runf/%s.tsimf > runf/%s.tsimf',     vhe_model, ident);  string= strcat(string,strin2);

#  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); exit; end

  string= sprintf('tsim -no_addarg -overwrite runf/%s.tsimf > qout/tsim_%s.out', ...
                                                  ident,      ident);
  [istat, output]= system(string); 
  if (istat ~= 0)
   printf(' ***ERROR*** In wear_loop.m\n');
   printf('             Exit status from program TSIM not equal 0 (zero)\n');
   printf('             Please open and read the error message in file:\n');
   printf('             qout/tsim_%s.out\n', ident);
   exit
  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); exit; end
  [istat, output]= system('del -q mplotr/wear_111r.print'); if (istat ~= 0); error(output); exit; end
  [istat, output]= system('del -q mplotr/wear_112l.print'); if (istat ~= 0); error(output); exit; end
  [istat, output]= system('del -q mplotr/wear_112r.print'); if (istat ~= 0); error(output); exit; end
  [istat, output]= system('del -q mplotr/wear_121l.print'); if (istat ~= 0); error(output); exit; end
  [istat, output]= system('del -q mplotr/wear_121r.print'); if (istat ~= 0); error(output); exit; end
  [istat, output]= system('del -q mplotr/wear_122l.print'); if (istat ~= 0); error(output); exit; end
  [istat, output]= system('del -q mplotr/wear_122r.print'); if (istat ~= 0); error(output); exit; 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)
   printf(' ***ERROR*** In wear_loop.m\n');
   printf('             Exit status from program MPLOT not equal 0 (zero)\n');
   printf('             Please open and read the error message in file:\n');
   printf('             qout/mplot_%s.out\n', ident);
   exit
  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
