{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 nmaxspeeds[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');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; }