1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:

//==========================================================================*
// HHR-Track: Korrekturen für Quadrat der maximal möglichen Geschwindigkeit
//--------------------------------------------------------------------------*
function TDriver.CalcHHRTrack(Sec: TSection; Value: Tdble): Tdble;
begin
  Result := Value;
  if (Sec.ID > 430) and (Sec.ID < 550) then // In beiden Fällen aktiviert
    Result := 1250 * Sec.SpeedSqrFactor
  else if (Sec.ID > 1500) and (Sec.ID < 1550) then //<<< Spitzkehre
    Result := 400 * Sec.SpeedSqrFactor;  // nur ohne Fkt. unten aktiviert!
end;
//==========================================================================*					

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:

//==========================================================================*
// Spitzkehre nach der Parabolika
//--------------------------------------------------------------------------*
procedure Spitzkehre
  (TD: TTrackDescription; SecS1, SecS2, SecE1, SecE2: TSection);
var
  ...
begin
  // Tangente Anfahrtsgerade mit Offset von 7 m nach außen bestimmen
  P1 := Add(SecS1.Mid2D,Mult(-7,SecS1.Tor2D));
  P2 := Add(SecS2.Mid2D,Mult(-7,SecS2.Tor2D));
  TS := TStraight.Create                         // Tangente aus zwei Punkten
    (P1.X,P1.Y                                   // der Anfahrtsgeraden
    ,P2.X-P1.X                                   // berechnen
    ,P2.Y-P1.Y);

  // Tangente Abfahrtsgerade mit Offset von 7 m nach außen bestimmen
  P1 := Add(SecE1.Mid2D,Mult(-7,SecE1.Tor2D));
  P2 := Add(SecE2.Mid2D,Mult(-7,SecE2.Tor2D));
  TE := TStraight.Create                         // Tangente aus zwei Punkten
    (P1.X,P1.Y                                   // der Abfahrtsgeraden
    ,P2.X-P1.X                                   // berechnen
    ,P2.Y-P1.Y);

  // Tangenten-Schnittpunkt bestimmen
  T := TS.Intersect(TE);

  // Tangenten-Schnittwinkel bestimmen
  Beta := ArcCos(TS.Angle(TE));

  // Zwangspunkt bestimmen (Scheitel der Kurve)
  Q := MinDist(TD,T,SecS2,SecE1,Ind);

  // Abstand des Zwangspunkts vom Tangentenschnittpunkt bestimmen
  TQ := Dist(T,Q);

  // Abstand des Zwangspunkts von Anfahrtstangente bestimmen
  Y := TS.Dist(Q);

  // Station des Zwangspunkts auf Anfahrtstangente bestimmen
  TAs := Sqrt(TQ * TQ - Y * Y);

  // Abstand des Punktes A auf der Anfahrtstangente von T bestimmen
  // Bei A beginnt der Kreisbogen
  Tgbh := tan(Beta / 2.0);
  YTgbh := Y * Tgbh;
  X := YTgbh + Sqrt(YTgbh*YTgbh +  2 * TAs * YTgbh - Y*Y);
  TA := TAs + X;

  // Radius bestimmen
  R := TA * tgbh;

  // Kontrolle:
  X2 := R*R - (R-Y)*(R-Y);
  if Abs(X2 - (X*X)) > 0.01 then
    WriteLn(Format('?: %.3f² = %.3f <> %.3f',[X, X*X, X2]));

  // Kreismittelpunkt schätzen
  DTQ := Normalize(Sub(Q,T));
  M := Add(T,Mult(R+TQ,DTQ));

  // Offsets der Sections bestimmen
  // 1. Kreisbogen vom Scheitel zurück bis etwa zum Punkt A
  Offset := FLT_MAX;
  Sec := TD[Ind];                                // Segment d. Zwangspunktes
  repeat
    LastOffset := Sec.Offset;
    Sec := TD[TD.Loop(Sec.ID - 1)];              // Vorangegangenes Segment
    Sec.Offset := NewOffset(Sec,M,R);
    if Offset > Sec.Offset then
      Offset := Sec.Offset;
  until (Sec = SecS2) or (Offset < Sec.Offset);
  Offset := Sec.Offset;

  // 2. Übergang zurück
  oP[0].X := 0.0;
  oP[0].Y := SecS1.Offset;
  oP[0].S := 0.0;
  oP[1].X := 60.0;
  oP[1].Y := -7;
  oP[1].S := 0.0;
  oP[2].X := 70.0;
  oP[2].Y := Offset;
  oP[2].S := 0.0;
  oP[3].X := 71.0;
  oP[3].Y := LastOffset;
  oP[3].S := 0.0;

  oOrdinate := TSplineOrdinate.Create;
  oOrdinate.Count := 4;
  oOrdinate.Points := Addr(oP);

  Delta := 70;
  repeat
    Sec := TD[TD.Loop(Sec.ID - 1)];              // Vorangegangenes Segment
    Delta := Delta - 1.0;
    Sec.Offset := oOrdinate.Evaluate(Delta);
  until Delta <= 0.0;
  oOrdinate.Free;

  // 3. Kreisbogen vom Scheitel vorwärts
  Offset := FLT_MAX;
  Sec := TD[Ind];                                // Segment d. Zwangspunktes
  repeat
    Sec := TD[TD.Loop(Sec.ID + 1)];              // Folgendes Segment
    Sec.Offset := NewOffset(Sec,M,R);
    if Offset > Sec.Offset then
      Offset := Sec.Offset;
  until (Sec = SecE1) or (Offset < Sec.Offset);
end;
//==========================================================================*					

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:

Original:
begin
  AllowedSpeed := GetAllowedSpeed(oCurrSec);
  //Gib Vollgas bis 1 m/s Unterschied
  if AllowedSpeed > SpeedX + 1.0 then
    Result := 1.0
  else
und Fälschung
  //Gib Vollgas bis 10% vom maximal erlaubten m/s Unterschied
  if AllowedSpeed/SpeedX >1.10 then
    Result := 1.0
  else					

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:

//==========================================================================*
// Ansteuerung eines Zielpunktes
//--------------------------------------------------------------------------*
function TDriver.GetTargetPoint(LookAheadDistance: Tdble): TV2D;
var
  ...
  Plan: Tdble;
begin
  if oPit.InPit then
  begin
    ...
  end
  else
  begin
    ...
  end;

  Plan := oMyOffset;                         // Aktuelle Querabeichung merken
  Offset := PlanOffset;                      // Querabweichung planen
  if Abs(Offset) > 0.01  then                // Wenn Planung nicht 0 ist
    Plan := Offset;                          //   die verwenden
  if oPit.InPit then                         // Weg in Boxengasse nehmen?
    Offset := Sign(Plan) * Dist(S,P);        //   Querabweichung berechnen
  Result := Add(S,Mult(-Offset,N));          // Zielpunkt berechnen
end;
//==========================================================================*					

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:

//==========================================================================*
// Aerodynamische Balance
//--------------------------------------------------------------------------*
procedure TFormMain.EditAngleChange(Sender: TObject);
var
  I: Integer;
  Hm: Double;            // ride height
  BodyLiftFront: Double; // Downforce at front
  BodyLiftRear: Double;  // Downforce at rear
  WingLiftFront: Double; // Downforce at front
  WingLiftRear: Double;  // Downforce at rear
  WingDragFront: Double; // Drag force at front
  WingDragRear: Double;  // Drag force at rear
  Ratio: Double;         // Ratio (Aerodynamische Balance)
  Line: string;
begin
  try
    MemoSpoiler.Lines.Clear;
    oCar.WingAngle[TEdit(Sender).Tag] := TEdit(Sender).Text;
    Hm := 0.0;
    // Mean weighted ride height at all wheels [m]
    for I := 0 to 3 do
      Hm := Hm + oCar.RideHeightVal[I];
    Hm := 1.5 * Hm;
    Hm := Hm * Hm;
    Hm := Hm * Hm;
    Hm := 2 * exp(-3.0 * Hm);

    BodyLiftFront := oCarType.Clift[FRONT] * HM;   // Down force at front
    BodyLiftRear := oCarType.Clift[REAR] * HM;     // Down force at rear

    WingLiftFront := 4.0 * 1.23                    // Down force at front
      * oCarType.WingArea[FRONT] * sin(oCar.WingAngleVal[FRONT]);
    WingLiftRear := 4.0 * 1.23                     // Down force at rear
      * oCarType.WingArea[REAR] * sin(oCar.WingAngleVal[REAR]);
    Ratio := (BodyLiftFront + WingLiftFront) /
      (BodyLiftFront + WingLiftFront + BodyLiftRear + WingLiftRear);

    Line := 'Resultierende aerodynamische Balance: ' + FloatToStr(Ratio);
      MemoSpoiler.Lines.Add(Line);

    WingDragFront := 1.23                          // Drag force at front
      * oCarType.WingArea[FRONT] * sin(oCar.WingAngleVal[FRONT]);
    WingDragRear := 1.23                           // Drag force at rear
      * oCarType.WingArea[REAR] * sin(oCar.WingAngleVal[REAR]);

    Line := 'Resultierender Luftwiderstandsbeiwert vorne: '
      + FloatToStr(WingDragFront);
    MemoSpoiler.Lines.Add(Line);
    Line := 'Resultierender Luftwiderstandsbeiwert hinten: '
      + FloatToStr(WingDragRear);
    MemoSpoiler.Lines.Add(Line);
  except
    Line := 'Ungültige Eingabe';
    MemoSpoiler.Lines.Add(Line);
  end;
end;
//==========================================================================*					

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:

    OldFactor := oCurrSec.SpeedSqrFactor;        // Aktueller Lernfaktor
    W := (Segment.Width * oWidthDiv - oBorder-CarWidth);
    Factor :=Abs(oActOffset)/W/oToleranz;
    {//Ausserhalb der Fahrbahn
    if Factor >EinEintel then
      begin
      Delta := (Factor-EinEintel)*oMMinus*OldFactor;
      N := 10;
      end
    else
    }
      begin
      if Factor < EinHalb then    // Je kleiner die Abweichung
        begin                     // umso größer die Erhöhung
        Delta := (Factor-EinHalb)*OldFactor *oMPlus;
        Writeln(oCurrSecID:5,'Schneller ',-delta:10:7,ABS(oActOffset):13:7,oActOffset:13:7,W:13:7);
        N := 6;
        end
      else //if oActOffset > oToleranz then          // bzw. die Erniedrigung
        begin                                        // Grenze ist 1,25 m
        Delta := (Factor-EinHalb)*OldFactor *oMminus;
        Writeln(oCurrSecID:5,'langsamer ',-delta:10:7,ABS(oActOffset):13:7,oActOffset:13:7,W:13:7);
        N := 12;
        end;
     end;
//...					

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:

  // Quer-Abweichung vom geplanten Verlauf feststellen
  oActOffset := oCurrSec.GetOffset(Pos); // Aktueller Wert, nur in Kurven <> 0!!!
  if oActOffset = 0.0 then               // Wenn ich auf einer Geraden bin, 
  begin                                  //   Dann nehme ich den Wert der letzten
    oActOffset := oLastOffset;           //   Kurve und lasse den langsam ausklingen
    oLastOffset := oLastOffset * 0.98;   //   um zu wissen ob ich nach außen
  end                                    //   abweiche oder nach innen
  else                                   // ansonsten
  begin                                  //   merke ich mir den aktuellen Wert
    oLastOffset := oActOffset;           //   der Kurve für die folgende Gerade
  end;					

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:

  Factor :=2.0*Abs(oMyOffset)/(Segment.Width*oWidthDiv{0.0..1.0}-CarWidth-oBorder)/oToleranz;
 //Sind sehr nah an der Mitte dann rauf
  Delta :=0;
  if Factor < EinHalb then    // Je kleiner die Abweichung
    begin                                        // umso größer die Erhöhung
    Delta := (Factor-Einhalb)*oMPlus;
    N := 3;
    end
  else
     begin
     if Factor > 1 then
        Delta := 2.0*oMMinus
     else 
       if Factor > 0.7 then
         //Nur Kurvenaußenseite wird gedrosselt 
         IF oMyOffset*oCurrSec.Rad<0.0 then
              Delta := (Factor-Einhalb)*oMMinus;
          //Writeln(oCurrSecID,' ',Factor:4:2,' zu weit: ',Delta:10:7);
          N := 4;
          end;
  IF Delta <> 0.0 then
    begin
    Delta := OldFactor*Delta/LookAhead;// *1/N weggelassen
    for I := 0 to Round(N * LookAhead) do        // Lernfaktor auf den letzen
      begin                                      // X m anpassen
      Sec := oTrackDesc[oTrackDesc.Loop(oCurrSecID - I)]; 
      tm := (Sec.SpeedSqrFactor - Delta);
      if (tm < 0.3) then
        tm := 0.3
      else
        If (tm>3.0) then
          tm := 3.0;
      Sec.SpeedSqrFactor := tm;
      end;
  end;
end;					

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:

  Fm  := mass*g;// Gewichtskraft

  for I := 0 to oCount - 1 do
  begin
    Sec := Section[I];
    Mu := Sec.oSeg.Surface.KFriction * TireMu * MuFactor;
    R := ABS(Sec.Rad);

    sinA :=sin(sec.okAlpha);
    cosA :=cos(sec.okAlpha);
    sinG :=sin(sec.okGamma);
    cosG :=cos(sec.okGamma);

    //Für kleine Winkel
    Fmq := -sinA*Fm; //Hangabtrieb quer
    Fml := sinG*Fm; //Hangabtrieb längs
    Fms := sqrt(sqr(Fm)-sqr(Fmq)-sqr(Fml)); //Der Rest senkrecht
    {
    Ff  := V*V*Abs(mass/R);   //Fliehkraft Krümmungskreis in Querrichtung
    Ffs := sinA*Ff;      //senkrecht zur Fahrbahn
    Ffq := cosA*Ff;      //quer zur Fahrbahn die eigentliche Wirkrichtung
    Ffl := 0.0;          //in Längsrichtung

    Fas := CA*V*V;         //Aerodynamische Anpressung genau senkrecht zur Fahrbahn

    Fs := Fms+Ffs+Fas;
    Fq := Fmq+FFq;
    //Fq==mue*Fs-> Grenzbedinung fürs Rutschen
    //mu*Fms-Fmq=Ffq-mu*Ffs
    //mu*...=V2*Abs(mass/R)*(cos(Sec.okAlpha)-mu*sin(Sec.okAlpha))
    //V2=oSSq=
    }
     
    IF R > 2000 then    //Sehr grosser Radius hat keinen Einfluß
      Sec.oSSq := MAX_SSq
    else
      begin
      //Neigung Strasse passend zur Kurve
      If Sec.Rad<=0.0 then
        V0 :=(cosA/mu+sinA)*mass/R-CA
      else
        V0 :=(cosA/mu-sinA)*mass/R-CA;

      IF V0 < 0 then //In überhöhten Kurven kann das passieren
        Sec.oSSQ := MAX_SSq
      else
        Sec.oSSq:= Fs/V0;
      end;
  end;					

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
154:
155:
156:
157:
158:
159:
160:
161:
162:
163:
164:
165:
166:
167:
168:
169:
170:
171:
172:
173:
174:
175:
176:
177:
178:
179:
180:
181:
182:
183:
184:
185:
186:
187:
188:
189:
190:
191:
192:
193:
194:
195:
196:
197:
198:
199:
200:
201:
202:
203:
204:
205:
206:
207:
208:
209:
210:
211:
212:
213:
214:
215:
216:
217:

program OptBremsWeg;
uses
  debug;
type
  Tdble = extended;
const
  g : tDble = 9.81;
  s_max = 10000.0;
var
  A,
  mue,
  alpha,
  sinA,
  cosA,
  gamma,
  sinG,
  cosG,
  yR,
  zR,

  mass,

  mueFahrbahn,
  mueRoll,
  mueSumm ,
  Luftdichte,
  Cw,
  ca,
  V0,
  V,
  V2,
  Ve,
  Vw,
  E,
  dE,


  F,
  Fl,
  Fq,
  Fs,
  Fb,
  Fmax,
  Flges,

  Ff,
  Ffq,
  Ffl,
  Ffs,

  Fm,
  Fmq,
  Fml,
  Fms,

  Fa,
  Faq,
  Fal,
  Fas,

  Fk,
  Fkq,
  Fkl,
  Fks,

  h,
  s,ds :Tdble;

function DeltaE(V2,ds:Tdble):Tdble;
begin
  //sekrechte Kräfte auf die geneigte Fläche umrechnen
  sinA := sin(Alpha);
  cosA := cos(Alpha);
  sinG := sin(Gamma);
  cosG := cos(Gamma);

  Fm  := mass*g;
  Fmq := cosA*sinG*Fm; //Hangabtrieb quer
  Fml := sinA*cosG*Fm; //Hangabtrieb längs
  Fms := sqrt(sqr(Fm)-sqr(Fmq)-sqr(Fml)); //Der Rest senkrecht

  Ff  := mass*V2/yR;   //Fliehkraft Krümmungskreis in Querrichtung
  Ffs := sinA*Ff;      //senkrecht zur Fahrbahn
  Ffq := cosA*Ff;      //quer zur Fahrbahn die eigentliche Wirkrichtung
  Ffl := 0.0;          //in Längsrichtung

  Fk  := mass*V2/zR;   //Kuppen und Täler,Krümmumgkreis in Längsrichtung

  Fks := cosG*Fk;      //senkrecht zur Fahrbahn die eigentliche Wirkrichtung
  Fkq := sinG*Fk;      //quer zur Fahrbahn
  Fkl := 0.0;          //in Längsrichtung

  Fas := ca*V2; //Aerodynamisches anpressen;// Wirkt genau senkrecht auf die Fahrbahn
  //Faq         //Aerodynamische Querkraft Seitenwind
  Fal := cw*V2; //Aerodynamisches Bremsen Luftwiderstand

  Fs := Ffs+Fms+Fks+Fas;  //Kraft senkrecht auf Fahrbahn
  Fq := Ffq+Fmq+Fkq;      //Kraft quer zur Fahrbahn
  Fl := Ffl+Fml+    Fal;  //Kraft in Fahrrichtung

  //Auf Fahrbahn durch die Reifen übertragbare Gesamtkraft
  IF Fs< 0.0 then
    Begin
    writeln(S,' Wir heben ab ',V);
    halt(-1);
    end;

  Fb := mueSumm*Fs;
  //Bremskraft des Wagens in Wagenrichtung ,Die Querkraft muss auch durch
  //die Reifen übertragen werden.
  Fmax := sqr(Fb)-sqr(Fq);
  //Rutschen wir schon, dann war es das
  IF Fmax > 0 then
    Fmax := sqrt(Fmax)
  else
    Begin
    Fmax := 0.0;
    writeln(S,' Wir rutschen ',V);//Aber Längskräfte bremsen oder beschleunigen noch
    halt(-1);
    end;

  Flges  := Fmax+Fl;
  {
  writeln('rel Beschleunigung : ',Flges/(mass*g));
  Writeln('Fg      : ',Fm);
  Writeln('Fflieh  : ',Ff);
  Writeln('FLoop   : ',Fk);
  Writeln('FAeroCa : ',Fas);
  Writeln('FAeroCw : ',Fal);
  }
  DeltaE := -Flges*ds;
end;

function Heun(V,h,ds:tdble):tdble;
//ausreichend genau, wenn sich nicht sich zu viele Parameter gleichzeitig
//während eines Schrittes aendern
var
  E0,dE1,dE2,
  h1,V2: tDble;
Begin
  V2 := sqr(v);
  E0 := (g*h +0.5*sqr(V))*mass;
  h := h+ds*sinA;//neue geodätische Höhe man kann auch Startziel als 0 nehmen
  dE1 := DeltaE(V2,ds);
  V2 := ((E0+dE1)/mass-g*h)*2.0;
  dE2 := DeltaE(V2,ds);
  V2 := ((E0+0.5*(dE1+dE2))/mass-g*h)*2.0;
  If V2 > 0 then
    V := sqrt(V2)
  else
    V:= 0;//Rollt eigentlich rückwärts, kann auch ein überschiessen sein
//  Writeln(s:8:4,E0:15:6,dE1:15:6,dE2:15:6,V:15:6);
  Heun := V;
end;

begin
  V0 := 30;
  A := 2.0;
  alpha :=-10.0*pi/180;   //Längsneigung
  gamma := 0*pi/180;   //Querneigung
  sinA := sin(Alpha);
  cosA := cos(Alpha);
  sinG := sin(Gamma);
  cosG := cos(Gamma);
  yR := 100;//Enge Kurve
  zR:= 1e38; //+Innenlooping(bergab->bergauf),bergauf -> bergab
  mueFahrbahn := 1.0;
  mueRoll:= 0.00;
  Luftdichte := 1.2;
  Cw := 0.5*Luftdichte*0.41*A;
  ca := 0.5*Luftdichte*2.3*A;
  mass := 1050.0;

  mueSumm := mueFahrbahn+mueRoll;


  ds := 8.0;
repeat;
  s := 0.0;
  h := 1000;
  V := V0;
  v2:=sqr(v);
  E := (0.5*V2+g*h)*mass;
  repeat
    h := h+ds*sinA;
    E := E+DeltaE(V2,ds);
    V2 := (E/mass-g*h)*2.0;
    If V2 < 0 then
      break;
    V := sqrt(V2);
    s := s+ds;
    If Abs(s-48) < ds then
      writeln(s:6:3,V:10:6,ds:10:6);
  until s>s_max;//Notbremse
  ds := ds*0.5;
until ds <1/8;
  writeln;


  ds := 8.0;
repeat
  s := 0.0;
  h := 1000;
  V := V0;
  E := (0.5*sqr(V)+g*h)*mass;

  repeat
    V:=Heun(V,h,ds);
    s := s+ ds;
    If Abs(s-48) < ds then
      writeln(s:6:3,V:10:6,ds:10:6);
  until (V <= 0.0) or (s>s_max);
  writeln(s:10:6,v:15:6,Heun(V,h,ds):15:6);
  ds := ds*0.5;
until ds <1/8;

end.					

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:

(Euler-Cauchy)
S      V          ds
48.000 22.980188  8.000000
48.000 22.610156  4.000000
48.000 22.416512  2.000000
48.000 22.317575  1.000000
48.000 22.267585  0.500000
48.000 22.242462  0.250000
48.000 22.229868  0.125000

(Heun)
48.000 22.273808  8.000000
48.000 22.232475  4.000000
48.000 22.221205  2.000000
48.000 22.218260  1.000000
48.000 22.217507  0.500000
48.000 22.217316  0.250000
48.000 22.217269  0.125000