{Use Pascal Pretty Printer, (C)2001 by Lavrov Dmyry} unit runge; interface type real=double; const maxlongint=\$7fFFffFF; const MaxRealArraySize=maxlongint div sizeof(real); type TrealArray=array[1..MaxRealArraySize]of real;pRealArray=^TrealArray; type Tuserfunc=procedure(t:real;n:longint;v,d:pRealArray); procedure initRK4(_f:TuserFunc;_n:longint;_dt:real;_p:pointer{pointer to array of real}); procedure StepRK4; var t:real; Phase:longint; dt:real; implementation var p:PrealArray; tmp1,tmp2,tmp3:PrealArray; f:TuserFunc; nEQN:longint; procedure initRK4(_f:TuserFunc;_n:longint;_dt:real;_p:pointer{pointer to array of real}); begin (**************disposing*******************) if tmp1<>nil then dispose(tmp1); if tmp2<>nil then dispose(tmp2); if tmp3<>nil then dispose(tmp3); nEQN:=_n; f:=_f; dt:=_dt; p:=_p; getmem(tmp1,sizeof(real)*nEQN); getmem(tmp2,sizeof(real)*nEQN); getmem(tmp3,sizeof(real)*nEQN); end; procedure StepRK4; var i:longint; begin if f=nil then exit; phase:=0; f(T,neqn,p,tmp1); {0} for i:=1 to neqn do tmp2^[i]:=p^[i]+dt*tmp1^[i]/2; inc(phase); f(T+dt/2,neqn,tmp2,tmp3);{1} for i:=1 to neqn do begin tmp1^[i]:=tmp1^[i]+2*tmp3^[i]; tmp2^[i]:=p^[i]+dt*tmp3^[i]/2; end; inc(phase); f(T+dt/2,neqn,tmp2,tmp3);{2} for i:=1 to neqn do begin tmp1^[i]:=tmp1^[i]+2*tmp3^[i]; tmp2^[i]:=p^[i]+dt*tmp3^[i]; end; inc(phase); f(T+dt,neqn,tmp2,tmp3);{3} for i:=1 to neqn do begin p^[i]:=p^[i]+dt*(tmp1^[i]+tmp3^[i])/6; end; t:=t+dt; end; begin tmp1:=nil; tmp2:=nil; tmp3:=nil; end.