{Use Pascal Pretty Printer, (C)2001 by Lavrov Dmyry}
{(C) 2002,2003 by Dmytry Lavrov.  http://dmytrylavrov.narod.ru}
{$mode delphi}
{define dump}
uses crt,GraphObj,ObjScr,p256,msmouse,
runge;
{------------------------------defs-----------------------------}
const KnutPoints=200;
KnutMass=1000;{KG}
KnutMaxForce=20000;{20 KN}
KnutK=100000{1/10 with force=10000 };
KnutLength=1e6;
pointLength=KnutLength/KnutPoints;
PointMass=KnutMass/KnutPoints;
CapsuleMass=pointmass;
TimeStep:real=0.1;
engFmaxY=4000;
engFmaxX=0{-3000};
engT2X=100;{time when engine has 50% of force}
engT2Y=100;{time when engine has 50% of force}
diss=0.0{01};{"dissipation"-very bad model of them}
{checking options}
mincosangle=0.9;


{show params}
CanBeSplited=true;
M=0.0002;{Марштаб}
speedshowM=0.01;
pointstoskip=0;
Calcs{BetweenRenders}=10;
FracParam=0.01;


{---------------------------------------------------------------}
{$ifdef dump}Var Dumpdata:array[1..100,1..100] of real;
DumpPhase:longint;{$endif}




Var SpecScreen:Tgraph;
var i,o:text;
type pnt=record
  x,y{,sx,sy}:real;
end;

var started:boolean;
Splited:boolean;

function EngineForce(t:real):pnt;
begin;
result.x:=0;
if started and(not splited) then
result.y:=(1-1/(1+sqr(t/engT2Y)))*engFmaxY
else result.y:=0;
end;
type TKnutpoint=record
  u,p:pnt;
end;
type Tpoints=record
data:array[1..KnutPoints]of TKnutpoint;
engine_energy:real;
end;
ppoints=^tpoints;
var CurPoints:Tpoints;

function sqrp(a:pnt):real;
begin
  sqrp:=sqr(a.x)+sqr(a.y);
end;
function Energy:real;
var n:longint;
r,d2:real;
lastP:pnt;
begin
  r:=0;
  lastP:=CurPoints.data[1].p;
  for n:=1 to KnutPoints do with curpoints.data[n] do begin
    r:=r+sqrp(u)/2;
    d2:=sqr(p.x-LastP.x)+sqr(p.y-LastP.y);
    if d2>pointLength*pointLength then r:=r+KnutK*sqr(sqrt(d2)-pointLength)/2;
    lastP:=curpoints.data[n].p;
  end;
  Energy:=r;
end;

var accel:array[1..KnutPoints]of pnt;

var splitedPts:array[1..KnutPoints] of boolean;
var MaxSpeeds:array[1..KnutPoints] of real;
var phase:boolean;
procedure MidSpeed;
var i:longint;
begin
if diss=0 then exit;
i:=ord(phase)+1;{1 or 2}
repeat
{bad model of dissipation}
curpoints.data[i].u.x:=(curpoints.data[i].u.x+curpoints.data[i+1].u.x*diss)/(1+diss);
curpoints.data[i].u.y:=(curpoints.data[i].u.y+curpoints.data[i+1].u.y*diss)/(1+diss);
curpoints.data[i+1].u.x:=(curpoints.data[i+1].u.x+curpoints.data[i].u.x*diss)/(1+diss);
curpoints.data[i+1].u.y:=(curpoints.data[i+1].u.y+curpoints.data[i].u.y*diss)/(1+diss);
i:=i+2;
until i>=KnutPoints;
phase:=not(phase);
end;
procedure Func(t:real;n:longint;v,d:pRealArray);
var i:longint;
F,dp:pnt;
d2,dl,fk:real;
begin
    F:=engineforce(t);
    ppoints(d)^.engine_energy:=F.x*ppoints(v)^.data[1].u.x+F.y*ppoints(v)^.data[1].u.y;
    ppoints(d)^.data[1].u.x:=F.x;ppoints(d)^.data[1].u.y:=F.y;{temporary,it's force!}
    ppoints(d)^.data[2].u.x:=0;
    ppoints(d)^.data[2].u.y:=0;
    for i:=1 to KnutPoints-1 do with ppoints(v)^.data[i] do begin
      dp.x:=p.x-ppoints(v)^.data[i+1].p.x;
      dp.y:=p.y-ppoints(v)^.data[i+1].p.y;
      d2:=sqr(dp.x)+sqr(dp.y);
      if not((splitedpts[i]or splitedpts[i+1])and CanBeSplited){and(d2>sqr(pointlength))} then
      begin
        dl:=sqrt(d2);
        fk:=(dl/pointLength-1)*knutK;
        if (fk>KnutMaxForce)and(runge.phase=0) then begin Splited:=true;SplitedPts[i]:=true; end;
        begin
          F.x:=(fk*dp.x/dl);{force of Pn and Pn+1}
          F.y:=(fk*dp.y/dl);
          ppoints(d)^.data[i+1].u.x:=F.x;{temporary it's force}
          ppoints(d)^.data[i+1].u.y:=F.y;
          ppoints(d)^.data[i].u.x:=(ppoints(d)^.data[i].u.x-F.x)/pointmass;{now it's acceleration}
          ppoints(d)^.data[i].u.y:=(ppoints(d)^.data[i].u.y-F.y)/pointmass;
        end;
      end else begin
      ppoints(d)^.data[i+1].u.x:=0;
      ppoints(d)^.data[i+1].u.y:=0;
      ppoints(d)^.data[i].u.x:=(ppoints(d)^.data[i].u.x)/pointmass;{now it's acceleration}
      ppoints(d)^.data[i].u.y:=(ppoints(d)^.data[i].u.y)/pointmass;
      end;
     ppoints(d)^.data[i].p.x:=ppoints(v)^.data[i].u.x;
     ppoints(d)^.data[i].p.y:=ppoints(v)^.data[i].u.y;
    end;{for    Acceleration cycle}
     inc(i);
     ppoints(d)^.data[i].u.x:=(ppoints(d)^.data[i].u.x)/capsulemass;
     ppoints(d)^.data[i].u.y:=(ppoints(d)^.data[i].u.y)/capsulemass;
     ppoints(d)^.data[i].p.x:=ppoints(v)^.data[i].u.x;
     ppoints(d)^.data[i].p.y:=ppoints(v)^.data[i].u.y;
end;
procedure recalc;
begin
  stepRK4;
end;
procedure initcalc;
begin
initRK4(@func,sizeof(Tpoints) div sizeof(real),timestep,@curpoints);
t:=0;
end;

function pnt2st(const p:pnt):string;
var tmp:string;
begin
str(p.x:3:3,tmp);
result:='[x='+tmp;
str(p.y:3:3,tmp);
result:=result+',y='+tmp+']';
end;
{$ifdef dump}
var j:longint;
procedure Dump;
var n:longint;
s:real;
begin
if DumpPhase>4 then begin
j:=j+1;
if j>100 then exit;
for n:=1 to 100 do with CurPoints.data[n*2] do begin
s:=sqrt(u.x*u.x+u.y*u.y);
dumpdata[n,j]:=s;
end;
DumpPhase:=0;
end;
inc(DumpPhase);
end;
procedure SaveDump;
var f:text;
n,m:longint;
begin
assign(f,'dump.txt');
Rewrite(f);
Writeln(f,'100x100');
for m:=1 to 100 do for n:=1 to 100 do writeln(f,dumpdata[n,m]);
close(f);
end;
{$endif}
var MaxSpeed, CAngle:real;
var cx,cy:longint;
procedure Render;
var n:longint;
iy:longint;
speed,force:real;
speedstep:real;
{angle checking}
var dx,dy,dl,lastdx,lastdy:real;
overangle:boolean;
{}
begin
  overangle:=false;
  with specscreen do begin
    clear(0);
    cur.x:=0;cur.y:=0;

    {for n:=1 to KnutPoints do}
    n:=1;
    dx:=0;
    dy:=0;
    while n<=KnutPoints do
    begin
      if SplitedPts[n] then curcolor:=colors[LightRed] else curcolor:=colors[15];
      {angle checking}
      if n<knutpoints then begin
      lastdx:=dx;
      lastdy:=dy;
      dx:=(curpoints.data[n+1].p.x-curpoints.data[n].p.x);
      dy:=(curpoints.data[n+1].p.y-curpoints.data[n].p.y);
      dl:=sqrt(dx*dx+dy*dy);
      if splitedPts[n]or splitedpts[n+1] then force:=0 else Force:=(dl/PointLength-1)*KnutK;
      dx:=dx/dl;
      dy:=dy/dl;
      cangle:=(dx*lastdx)+(dy*lastdy);
      if  cangle<=mincosangle then line(n div(pointstoskip+1),max.y div 2 +2,n div(pointstoskip+1),max.y div 2 +4);
      end;

      speed:=sqrt(sqr(curpoints.data[n].u.x)+sqr(curpoints.data[n].u.y));
      line(n div(pointstoskip+1),max.y div 2,n div(pointstoskip+1),round((max.y div 2)-speed*speedshowM));
      if speed>maxspeeds[n] then begin maxspeeds[n]:=speed; end;

      putpixel(n div(pointstoskip+1),round( (max.y div 2)-maxspeeds[n]*speedshowM),colors[10]);
      if force<0 then begin force:=-force;curcolor:=colors[LightBlue];end;
      {force}line(n div(pointstoskip+1),max.y-round((Force/KnutMaxForce)*Max.Y/2),n div(pointstoskip+1),max.y);

      if (speed>MaxSpeed)and(n>5) then begin MaxSpeed:=speed;curcolor:=10; end;
      PutPixel(round(curpoints.data[n].p.x*m)+cx,cy+round(curpoints.data[n].p.y*m));
    inc(n,pointstoskip+1);
    end;
    if Splited then writeln(o,'Нить порвана!(сила>',KnutMaxForce/1000:3:3,'KN)');
    Writeln(o,'time=',t:3:3,' ef=',pnt2st(engineforce(t)),' Max speed=',MaxSpeed/1000:3:3,' km/s ');
    n:=0;
    curcolor:=colors[7];
    SetWriteMode(xorPut);
    {speedshowM*speedstep>16}
    speedstep:=exp((round( ln(4*16/speedshowM)/ln(10) ))*ln(10))/4;
    repeat
    iy:=round((max.y div 2)-n*speedstep*speedshowM);
    if iy<20 then break;
    line(0,iy,(KnutPoints)div(pointstoskip+1)+1,iy);
    moveto((KnutPoints)div(pointstoskip+1)+2,iy-chrheight(#255) div 2);
    Write(o,n*speedstep/1000:3:3);
    n:=n+1;
    until iy<20;
    line(KnutPoints div 2,0,KnutPoints div 2,Max.Y div 2);
    curcolor:=colors[15];
    SetWriteMode(normalput);

        screen.PutImage(specscreen,0,0);
    screen.MoveTo(0,0);
  end;
end;
var rk:char;n:longint;  x,y:real;
var depth:longint;
procedure PrepareLine(s,e:longint);
var n:longint;
const nvaves=1;
begin
for n:=s to e do with curpoints.data[n].p do begin
y:=-sin((n/knutpoints)*2*nvaves*pi)*knutlength/80/nvaves;
end;
end;
{var m:longint;my:real;
begin
inc(depth);
if depth>80 then screen.outtext('????????');
if s<e-1 then begin
m:=(s+e)div 2;
my:=-(curpoints.data[s].p.y+curpoints.data[e].p.y)/2;
my:=my+(random-0.5)*pointlength*fracparam*(e-s);
curpoints.data[m].p.y:=my;
prepareLine(s,m);
prepareLine(m,e);
end;
dec(depth);
end;}


var imagename:string;
begin
{$ifdef dump}j:=0;{$endif}
  randomize;
  MaxSpeed:=0;
  InitScreen(0,0,0);screen.CurColor:=colors[15];screen.fontback:=1;
  specscreen.init;
  specscreen.AllocateMem(screen.max.x,screen.max.y);
  specscreen.AssignCrt(i);Reset(i);
  specscreen.AssignCrt(O);Rewrite(O);
  initcalc;
  cx:=screen.max.x div 2;
  cy:=screen.max.y div 2+10;

  curpoints.data[1].p.y:=0;
  curpoints.data[KnutPoints].p.y:=0;
  screen.outtext('line ');
  prepareline(1,KnutPoints);
  screen.outtext('ok');
{  x:=0;y:=0;
  for n:=1 to KnutPoints do with curpoints.data[n] do begin
    SplitedPts[n]:=false;
    x:=x+sqrt(sqr(pointlength)-sqr(p.y-y));
    p.x:=x;
    y:=p.y;
    u.x:=0;
    u.y:=0;
    {p.y:=sin((n/knutpoints)*5*pi)*knutlength/100;}
    MaxSpeeds[n]:=0;
    {u.y:=(abs((n-KnutPoints div 2)*10000/knutpoints)-2500)/8;}
  end;  }
  {  for n:=1 to KnutPoints div 2 do with curpoints.data[n] do begin
  u.x:=2000;
  u.y:=n-;
  end;}
  repeat
    for n:=1 to calcs do recalc;
    if not started then t:=0;
    {!!!!!!!!!!!}MidSpeed;{!!!!!!!!!!}
    render;

    {$ifdef dump}if started then Dump;{$endif}


    if keypressed then begin
      rk:=readkey;
      case rk of
        #0:begin rk:=readkey; end;
        #13,'s':started:=not started;
        't':begin screen.moveto(0,screen.max.y*3 div 4);write('timestep=');readln(TimeStep);runge.dt:=timestep;end;
        'c':begin Writeln;Write('Image name>');readln(imagename);if imagename<>'' then specscreen.writeimage(imagename+'.pcx');end;
      end;
    end;
  until rk=#27;
  SetTxMode;
  {$ifdef dump}
  SaveDump;
  {$endif}
end.






{------------------------------Voyager defs-----------------------------}
{const KnutPoints=200;
KnutMass=10000;{KG}
KnutMaxForce=200000{10000};{200 KN}
KnutK=1000000{1/10 with force=100000 };
KnutLength=1e6;
pointLength=KnutLength/KnutPoints;
PointMass=KnutMass/KnutPoints;
CapsuleMass=800;
TimeStep:real=0.1;
engFmaxY=80000;
engFmaxX=0{-3000};
engT2X=10;{time when engine has 50% of force}
engT2Y=10;{time when engine has 50% of force}
diss=0.0{01};{"dissipation"-very bad model of them}
{checking options}
mincosangle=0.9;
}