Помогите решить задачу, Написал программу для расчета температурной зависимости угольной части Опции

[Александр89]

Здравствуйте прошу помочь изменить код программы. Постановка задачи такова: Частица угля (шар) начинает прогреваться на границе действует лучистый и конвективный теплообмен. Решается методом сеток. По мере нагревания с температуры 400 градусов происходит реакция разложения. Проблема в том что программа вводит значения температуры от радиуса и степень разложения в конечный момент времени. А нужно что бы в каждый момент времени при выполнении условия 400 градусов считалась степень разложения и выводилась в файл( тоесть надо получить график изменение степени разложения от температуры) при этом на каждом промежутке степень разложения суммировалась. Вот код:

uses crt;

const mf=1000;
         sigma=5.669e-8;
         eps=1e-5;
         R=8.31;


type
     vector=array[1..mf] of real;
    var
 i, j, N  : integer;
 T, Tn,Tr: vector;
 alfa, beta,v: vector;
 ai, bi, ci, fi, d: real;
 lamda, ro, c, E, k0 : real;
 kapa, Te, eps1: real;
 h, tau, t_end, time : real;
 T0, R1, qxim : real;
 f, g : text;
begin
 clrscr;
 Writeln('   , N');
 Readln(N);
   Writeln('  , t_end');
 Readln(t_end);
 Writeln('  , R1');
 Readln(R1);
  Writeln(' , lamda');
 Readln(lamda);
 Writeln(' , ro');
 Readln(ro);
 Writeln('  , c');
 Readln©;
 Writeln('  , kapa');
 Readln(kapa);
 Writeln(' , Te');
 Readln(Te);
 Writeln('  , eps1');
 Readln(eps1);
 Writeln('  , T0');
 Readln(T0);
 Writeln('  , k0');
 Readln(k0);
 Writeln('  , qxim');
 Readln(qxim);
 Writeln('  , E');
 Readln(E);
 h:=R1/(N-1);
 tau:=t_end/1000.0;
 for i:= 1 to N do
      T[i]:=T0;
 time:=0;

 while time<t_end do
  begin
   time:=time+tau;
   for i:=1 to N do
   if T[i]<=373 then
   v[i]:=0
   else
    v[i]:=(exp(-k0*exp(-E/(R*T[i]))*tau));
 { for i:=1 to N do

  writeln(v[i]);            }

      begin
         for i:=1 to N do
     alfa[1]:=1.0;
     beta[1]:=0.0;
      for i:= 2 to N-1 do
       begin
      ai:=(lamda/(h*h))+(lamda/((i-1)*h*h));
      ci:=(lamda/(h*h))-lamda/((i-1)*h*h);
      bi:=ai+ci+ro*c/tau;
      fi:=(-ro*c*T[i]/tau)-qxim*v[i];
      alfa[i]:=ai/(bi-ci*alfa[i-1]);
      beta[i]:=(ci*beta[i-1]-fi)/(bi-ci*alfa[i-1]);
       end;
       repeat
       d:=T[N];
       T[N]:=(lamda*beta[N-1]+h*kapa*Te+eps1*sigma*h*((sqr(sqr(Te)))-(sqr(sqr(d)))))/(h*kapa+lamda*(1-alfa[N-1]));
       until abs(d-T[N])<=eps;
      for i:= N-1 downto 1 do
     T[i]:=alfa[i]*T[i+1]+beta[i];
   end;
     end;

   Assign(f,'res.txt');
 Rewrite(f);
 Writeln(f,' R1 = ',R1:6:4);
 Writeln(f,'  N = ',N);
  Writeln(f,'  lamda =  ',lamda:6:4);
 Writeln(f,'  ro = ',ro:6:4);
 Writeln(f,'   c  = ',c:6:4);
 Writeln(f,' T0 = ',T0:6:4);
 Writeln(f,'   kapa = ',kapa:6:4);
 Writeln(f,'  Te = ',Te:6:4);
 Writeln(f,' eps = ',eps:6:4);
 Writeln(f,' xh = ',h:6:4);
 Writeln(f,'  tau = ',tau:6:4);
 Writeln(f,'   t = ',t_end:6:4);
 close(f);
 Assign(g,'tempr.txt');
 Rewrite(g);
 for i:=1 to N do
    writeln(g,'  ',h*(i-1)8,'    ',T[i]:8:5, ' ', v[i]:10:10);
 close(g);
end.