Программа расчёта нагрева рубашки плазмотрона

#include <iostream>
#include <iomanip>
#include <vector>

int main()
{ const unsigned int LayersNumber=10;
  std::vector <float> v(LayersNumber); // Массив представляющий слои стенки
  const unsigned int TimeAccuracy=1000000;
  float S=.010, // Толщина
        T0=20, T1=20, a1=4e4, T2=3000, a2=1.5e3, TimeEnd=3;
  float HeatCapacity=3.44e4; // Объёмная теплоёмкость - проверить!
  float LayerHC=HeatCapacity*S/LayersNumber; // Теплоёмкость слоя
  float TermalConduction=401;
  float TCL=TermalConduction/S*LayersNumber; // Теплопроводность слоя
  float TimeStep=TimeEnd/TimeAccuracy;
  for (int i=0; i<LayersNumber; ++i) // Инициализация слоёв
  { v[i]=T0+(T1-T0)/LayersNumber*(i+.5); // Линейное распределение температуры
    std::cout << v[i] << " ";
  }
  std::cout << std::endl << "===========================" << std::endl;
  for (int t=0, f=1; t<TimeAccuracy; ++t) // Интегрирование теплоты по времени
  { for (int i=0; i<LayersNumber; ++i) // Проход по слоям
    { float Plus,  //Теплота поступающая от предыдущего слоя
            Minus; //Теплота отдаваемая следующему слою
      if (i==0) Plus=a2*(T2-v[i]);
      else Plus=TCL*(v[i-1]-v[i]);
      if (i<LayersNumber-1) Minus=TCL*(v[i]-v[i+1]);
      else Minus=a1*(v[i]-T1);
      v[i]+=(Plus-Minus)*TimeStep/LayerHC;
      // Корректировка грубого расчёта температуры нагрева
      // при большом шаге интегрирования по времени:
      if (i==0) if (v[i]>T2)     v[i]=T2;
                else ;
      else      if (v[i]>v[i-1]) v[i]=v[i-1];
                else ;
    }
    //if (t%(TimeAccuracy/10)==0) { for (int i=0; i<LayersNumber; ++i) std::cout << v[i] << " "; std::cout << std::endl; }
    if (t==f*10) f=t; // Изменение счётчика для увеличения периода вывода промежуточных расчётов
    if (t%f==0) // Промежуточный Вывод температуры слоя
    { std::cout << "/" << TimeEnd/TimeAccuracy*(t+1) << "/ ";
      for (int i=0; i<LayersNumber; ++i) std::cout << v[i] << " ";
      std::cout << std::endl;
    }
  }
  std::cout << "===========================" << std::endl;
  for (int i=0; i<LayersNumber; ++i) // Вывод температуры
  { //v[i]=T0+(T1-T0)/LayersNumber*(i+.5); // Линейное распределение температуры
    std::cout << v[i] << " ";
  }
  std::cout << std::endl << "===========================" << std::endl;

  //Оценочные проверки мощности и возможного нагрева, стационар
  std::cout << "Входной, проходной и  выходной потоки " << (T2-(v[0]+(v[0]-v[1])/2))*a2 <<
        " " << (v[0]-v[LayersNumber-1])*TermalConduction/(S*(LayersNumber-1)/LayersNumber) <<
        " " << (v[LayersNumber-1]-(v[LayersNumber-2]-v[LayersNumber-1])/2-T1)*a1 << std::endl;
  float A,B,C,D,E,x,y;
  A=T2;
  B=a2;
  C=TermalConduction/S;
  D=a1;
  E=T1;
  x=(A*B*(D+C)+C*E*D)/((B+C)*(D+C)-C*C);
  y=(E*D+C*x)/(D+C);
  std::cout << "Температуры поверхностей " << x << " " << y  << std::endl;
  std::cout << "Балланс " << (T2-x)*a2 << " " << (x-y)*TermalConduction/S << " " << (y-T1)*a1 << std::endl;
}