Program Plot_3D;

{ PLOT 3-D

  Ein Programm von Jens Gelhar, (c) Himpelsoft 1990.

  Unter Benutzung einiger Ideen (mein Programm kommt aber OHNE
  Quasi-Raytracing aus!) von Michael Aumer (KICKSTART 9/90, S.72)

  Dieses Programm plottet dreidimensionale Funktionen (z = f(x,y)).
  Dabei werden aber nicht irgendwelche häßliche Drahtgitter erzeugt
  (womöglich gar ohne Hidden Lines), sondern der Graph wird als
  plastischer Körper einschließlich Schatten dargestellt.

  Der Funktionsterm und einige Parameter müssen im Quelltext
  angegeben werden.

  Viel Spaß beim Ausprobieren!

}

{$opt q }

Uses Intuition, Graphics;

{
Durch Definition eines Bezeichners "test" wird die "Preview"-
Version des Plotters definiert, d. h. die Grafik wird erheblich
schneller, aber sehr grob erstellt:

Const test=1;

}

Const
  ScrH = 256;
  WinH = ScrH;
  ScrB = 640;
  WinB = ScrB;
  MinY = 80;
  MaxY = WinH-8;

  Schattenlänge = 4.5;   { größerer Wert -> längere Schatten }
  Schattenhärte = 20;    { größerer Wert -> weichere Schatten }

  {$if def test }
    Step = 15;
  {$endif }

Var scr: ^Screen;
    win: ^Window;
    ras: ^RastPort;
    vp : ^ViewPort;
    Msg : Ptr;
    xi, zi : integer;
    yi : real;
    YStep : real;
    x, y, z : real;
    i, pen : integer;
    xd, yd, maxYd : real;
    xu, yu: integer;
    v0, v1, v2 : real;
    cosA : real;
    lastZ : real;
    Schatten, Schattenlinie : integer;
    Txt: String;
    tf: real;
    max : Array[0..WinB] of integer;
    plotIt : Boolean;
    SL2 : Real;
    xj : integer;


{ Der folgende Text wird über die fertige Grafik geschrieben: }

Const Funktion = 'z = cos(c*r) / (r^2+1)';


Function f(x, y: real): real;
  { Diese Function repräsentiert die dargestellte Funktion.
    Wertebereich von x und y: -1 bis +1, wenn man andere Ausschnitte
    darstellen will, sind die Werte innerhalb der Function
    entsprechend zu normieren.
    Der von "f" zurückgegebene Wert muß der gewünschten Höhe in
    Bildpunkten (!) entsprechen. Die meisten Funktionen wird man also
    mit einem Faktor (z. B. 50) multiplizieren müssen, sonst wird das
    Bild zu "flach".
  }
  Var r : real;
  Begin
    r := sqr(x) + sqr(y);
    f := 50 * (1+cos(20*sqrt(r))) / (20*r+0.5);
  End;


(******************* andere schöne Funktionen: ******************]

  { Kreis mit Mittelpunkt: }
  Var r : real;
  Begin f := 0;
    r := sqr(x) + sqr(y);
    If r < 0.1 Then
      f := 80 * cos(5*r*pi)
    Else
    If abs(r-0.9) < 0.1 Then
      Begin
        r := r-0.9;
        f := 25 * cos(5*r*pi)
      End
    Else f := 0;
  End;

  { Überlagerung zweier gedämpfter Schwingungen: }
  Var r1, r2: Real;
  Begin
    r1 := 3*(sqr(x+0.5) + sqr(y+0.2));
    r2 := 3*(sqr(x-0.6) + sqr(y-0.5));
    f := 10 *(cos (8*r1) / (r1 + 0.05) + cos(8*r2) / (r2+0.05) );
  End;

  { "Chaotische" Wellen: }
  Begin
    f := 50 * ( sin (3*x+4*y) * abs(cos (3*y-2*x))) + 20 * (sqr(x) + 1-sqr(y) );
  End;

  { total "hohl": }
  Var r : real;
  Begin
    r := sqr(x) + sqr(y);
    If r < 0.9 Then f := -100*cos(r*pi/1.8)
             Else f := 0
  End;

[****************************************************************)


Function Tiefenfaktor (y: real): real;
  { perspektivische Verkleinerung, 0 < y < 1, normiert auf tf(0)=1 }
  Begin
    Tiefenfaktor := 2 / (y+2)
  End;


Begin
  {$opt b+ - Abbruch möglich }
  scr := Open_Screen(0, 0, ScrB, ScrH, 4, 0, 15, HIRES, 'Plot');
  win := Open_Window(0, 0, WinB, WinH, 1, MOUSEBUTTONS, BORDERLESS, Nil, scr,
                     0, 0, WinB, WinH);
  ras := win^.RPort;
  vp := ^scr^.ViewPort;

  For i:= 0 to 15 Do
    SetRGB4(vp, i, i, i, i);
  For i:= 0 to WinB do
    max[i] := MaxY;

  YStep := 2 / (MaxY-MinY);
  SL2 := Schattenlänge/2;

  yi := MaxY;
  While yi >= MinY Do
    Begin

      yu := Round(yi);

      y := (MaxY - yi) / (MaxY-MinY);       { [0..1] }
      y := 2*(exp(y)-1) / (exp(1)-1) - 1;   { [-1..+1] }

      tf := Tiefenfaktor( (MaxY-yi) / (MaxY-MinY) );

      Txt := IntStr(Round(50*y)+50) + '%';
      Move (ras, 1, 10);
      SetAPen (ras, 12);
      If GfxText (ras, Txt, Length(Txt)) Then;

      Schatten := -MaxInt;
      lastZ := f(-1-YStep, y);
      maxYD := 0;

      xi := 0;
      While xi < WinB Do
        Begin
          x := (2*xi - WinB) / WinB;
          xu := (WinB +Round( WinB * x * tf) )div 2;
          z := f(x,y);
          zi := yu - round( z*tf );

          plotIt := zi < max[xu];

          If plotIt Then
            Begin
              xd := (z - lastZ) {$if def test} / step {$endif} ;
              yd := z - f(x, y-YStep);

              { Vektorprodukt:  V := (yStep, 0, xd) X (0, yStep, yd) }
              v0 := -xd * YStep;
              v1 := -yd * YStep;
              v2 := sqr(YStep);

              { A = (-SL2, 0, 1),
                cos alpha = (A*V) / (abs(A) * abs(V)) }

              cosA := (v2 - SL2*v0) / (sqrt(1+sqr(SL2)) * sqrt(sqr(v0)+sqr(v1)+sqr(v2)));

              pen := round((14*16 div 2) * (cosA+1));

            End;

      Schattenlinie := xu + Round( (yu-zi) * Schattenlänge );
      If (Schattenlinie < Schatten) and plotIt Then
        Begin
          If Schatten-Schattenlinie < Schattenhärte Then
            pen := round(pen / (1+2*(Schatten-Schattenlinie)/Schattenhärte))
          Else
            pen := pen div 3
        End
      Else
        If Schattenlinie >= Schatten Then
          Schatten := Schattenlinie;

      If plotIt Then
        Begin

          {$if def test }
            pen := pen shr 4;
          {$else }
            If random(16) < pen and 15 Then
              pen := (pen shr 4) + 1
            Else
              pen := pen shr 4;
          {$endif }

            If yi = MaxY Then pen := 10;
          SetAPen (ras, pen);

          {$if def test }
            For i := 0 to Step-1 Do
              Begin
                xj := xu + round(tf*i);
                Move (ras, xj, zi+round((Step-i)/Step * (z-lastZ) * tf ));
                Draw (ras, xj, max[xj]-1);
              End;
            For i := 0 to Step-1 Do
              Begin
                xj := xu + round(tf*i);
                max[xj-1] := zi + round((Step-i)/Step * (z-lastZ) * tf )
              End;
          {$else }
            Move(ras, xu, zi);
            Draw(ras, xu, max[xu]-1);
            max[xu] := zi;
            If yd > maxYD Then maxYD := yd;
          {$endif }

        End;

          LastZ := z;

          xi := xi + {$if def test} step {$else} 1 {$endif} ;

        End;

     {$if def test}
       yi := yi - 1 - step div 2;
     {$else }
       If maxYd*tf > 3 Then             { bei starken Gefällen zum }
         yi := yi - 0.1 - 2/maxYd/tf    { Betrachter hin um weniger }
       Else                             { als eine Zeile springen }
         yi := yi - 1
     {$endif}

    End;

  Txt := Funktion;
  Move (ras, 1, 10);
  SetAPen (ras, 12);
  If GfxText (ras, Txt, Length(Txt)) Then;

  Repeat
    { alle bisherigen Events überlesen }
    msg := Get_Msg(win^.UserPort);
  Until msg = nil;

  Repeat
    { Programm mit linker Maustaste beenden }
    msg := Wait_Port(win^.UserPort);
  Until msg <> nil;

  { Fenster + Screen schließen sich autmatisch }

End.