#!/usr/bin/perl

#################                                       #########################
# Die Form der Schaufeln des Wasserrades vs. der Kurve des eintretenden Wassers #
# Skript zum ansteuern des Programms Plotmtv
# dies wird mit einer temporären Datei gefüttert
# Es werden default werte für die Radabmessungen vorgegeben
# Die Eingabe soll später menue geführt sein; zuerst einfach stdin
#
# Autor: Manfred Brücher mani.b@t-online.de
# 9.12.2001 
# sub prompt modified
# 20.1.2002
#
#### License: GPL ####
# 

use Math::Trig;

$g                = 9.807;
$r_Rad            = 2.5;
$r_Schaufel       = 1;
$r_innen_Schaufel = 0.085;
$grad_Schaufel    = 23;
$radiale_Tiefe    = 0.33;
$Schaufelzahl     = 54;
$v_Rad            = 1.57;     # Radumfangsgeschwindigkeit
$v_Wasser         = 3.15;       # Wassergeschwindigkeit   
$V_Wasser         = 0.3;     # Wassermenge
$b_Rinne          = 1.3;
$s_Rinne          = 0.01;    # der Luftspalt zwischen Rad und Rinne
$x_Rinne          = -0.0;    # Der negative horizontale Abstand der Rinne zum Radscheitel
$grad_Rinne       = 0;       # Die Richtung des Wassers auf der Rinne

$h_Wasser         = $V_Wasser/$b_Rinne/$v_Wasser; # Die dicke des eintretenden Strahls
$Radbreite        = 1.5;
$t_Rad            = 0.18; #-0.08;       # Der Zeitpunkt der Betrachtung der Drehung des Rades von x=0 an
#$start_drehung    = 0;       # Die Anfangsstellung des Rades in rad!!!! wird mit t_Rad berechnet
$t_Wasser	  = 0.00;     # Der Zeitpunkt der Betrachtung
$dt_Wasser        = 0.025; #0.32; 0.05;    # Wasserstrahllaenge; der delta Wert der Zeit
$t_schritt        = 0.020;    # die Zeitabstände in denen das Rad sich weiter dreht
$z_schritt        = 15;       # die Anzahl der Momentaufnahmen
$sichtbar         = 5;      # nur -5 - +5 Schaufeln berechnen

$i                = 1;       # Schleifenzähler nur zum debugen
$rad_schritt      = pi() / $Schaufelzahl / 1;  # in diesen Winkelschritten erfolgt die Berechnung (Genauigkeit)

$rad_S_teilung    = 2*pi()/$Schaufelzahl;
$rad_Schaufel     = deg2rad($grad_Schaufel);
$rad_Rinne        = deg2rad($grad_Rinne);

sub prompt
  { my ($einheit, $altwert, $wert);
    $einheit = $_[2] ne '' ? $_[2] : "m";
    print ("\n$_[0] \[ $_[1] $einheit \] \>");
    $altwert = $_[1];  
    if ($_[3] eq "1")   # 3. parameter 1 also auf eingabe warten
      { for (;;) 
         {$wert = <STDIN>;
	  chop($wert);
	  $wert =~ s/,/\./; # aus komma wird punkt
	  $wert =~ s/\s//g; # alle leerzeichen löschen
	  if ($wert =~ /\A-?\d*\.?\d*\z/)  # Zeilenanfang; 0-1 '-'; 0- dezimalstellen; 0-1 punkt; 0- dez.; Zeilenende
	   { if  ( $wert eq "") 
	       {$wert = $altwert;}
	     last; #return $wert; es kann nur ein return geben?!
	   }
	   print "\nNur Zahlen 0-9 dezimalpunkt \".\" und minus \"-\" sind erlaubt\nNochmal:> ";
	   next; 
	 }
     }else {$wert = $altwert;}
    #print "\n";
    return $wert;   # den Inhalt zurückgeben!
  }
  
sub nein
  { my ($w);
    print "@_ j/N > ";
    $w = <STDIN>;
    chop $w;
    $w = $w eq '' ?  "n" : $w;
    $w = $w eq 'N' ? 'n' : $w;
    if ($w eq "n"){return 1;};
   }


print ("\nDieses Skript versorgt das Programm plotmtv mit daten für das Wasserrad\n");
print ("zurzeit mit default Werten\n");

$mtv = "/tmp/wasserrad" . "$$.mtv";
print ("Die Ausgabe erfolgt in temporäre datei: $mtv\n");
open (MTV, "> $mtv")|| die ("Öffnen der tmp file nicht gelungen!");

print "Die aktuellen Parameter sind folgende:\n";
$m = 0;          # 0 oder ... nur ausgeben; 1 auf neue eingabe warten
do {
   print "\nDie Rad Abmessungen";
   $r_Rad         = prompt ("Der Radius       ", $r_Rad ,'' , $m);
   $radiale_Tiefe = prompt ("Radiale Tiefe    ", $radiale_Tiefe ,'' , $m);
   $Schaufelzahl  = prompt ("Anzahl Schaufeln ", $Schaufelzahl, "St.", $m);
   $Radbreite     = prompt ("Radbreite        ", $Radbreite, "m", $m);
   $b_Rinne       = prompt ("Breite der Rinne ", $b_Rinne, "m", $m);
   $m = 1;
   }
until  nein ("\nSollen die Radgrößen verändert werden?");
$m = 0;
do {   
   print "\nDie Schaufelform";
   $grad_Schaufel    = prompt ("Schaufelwinkel   ", $grad_Schaufel , "°" , $m);
   $r_Schaufel       = prompt ("Schaufelradius   ", $r_Schaufel    , "m" , $m);
   $r_innen_Schaufel = prompt ("Innerer Radius   ", $r_innen_Schaufel, "m", $m);
   $m = 1;
   }
until  nein ("\nSoll die Schaufelform verändert werden?");
$m = 0;
do {
    print "\nDie Wassermenge, Umfangsgeschwindigkeit ..., der Zeitpunkt der Betrachtung";
    $V_Wasser      = prompt ("Die Wassermenge    ", $V_Wasser, "m³/s" ,$m);
    $v_Rad         = prompt ("Die Umfangsgeschw. ", $v_Rad, "m/s", $m);
    $v_Wasser      = prompt ("Die Geschw. d. Wassers", $v_Wasser, "m/s", $m);
    $x_Rinne       = prompt ("Horiz. Abst. Radscheitel - Rinne ", $x_Rinne, "m", $m);
    $grad_Rinne    = prompt ("Die Neigung der Rinne geg. Horiz.", $grad_Rinne, "°", $m);
    $t_Rad         = prompt ("Zeitpunkt Rad      ", $t_Rad, "s" , $m);
    $m = 1;
    }
until nein ("\nWas ändern?");

print "\nRechne...\n";
$h_Wasser         = $V_Wasser/$b_Rinne/$v_Wasser; # Die dicke des eintretenden Strahls neu berechnen
$rad_S_teilung    = 2*pi()/$Schaufelzahl;
$rad_Schaufel     = deg2rad($grad_Schaufel);
$rad_Rinne        = deg2rad($grad_Rinne);

for ($z = 0; $z < $z_schritt; $z++)
{    $t_R      = $t_Rad    + $z * $t_schritt;
     $t_W      = $t_Wasser + $z * $t_schritt;
     
                                       #  boundary = off 
     print MTV "\$ DATA=CURVE2D\n\% toplabel = \"Wasserrad\"\n\% equalscale = on xmin = -0.7 xmax = 1.2 ymin = 1.6 ymax = 2.7\n";

     print MTV "\% linetype = 3 linewidth  = 0  linecolor = 1  markertype = 0  markercolor = 1\n";
     $bogen = 0; # den äußeren Umriss des Rades
     while ($bogen < 2 * pi() + $rad_schritt) {
	    $x = sin ($bogen) * $r_Rad;
	    $y = cos ($bogen) * $r_Rad;
	    print MTV "$x  $y\n";
	    $bogen = $bogen + $rad_schritt;
	    }


     print MTV "\n\n\%linetype = 1 linewidth = 1 linecolor = 1\n";
     $r = $r_Rad - $radiale_Tiefe; # den Radboden
     $bogen = 0;
     while ($bogen < 2 * pi() + $rad_schritt) {
	    $x = sin ($bogen) * $r;
	    $y = cos ($bogen) * $r;
	    print MTV "$x  $y\n";
	    $bogen = $bogen + $rad_schritt;
	    }

     $start_drehung = $v_Rad * $t_R  / $r_Rad;  # Der Winkel in radiant zum Zeitpunkt t_Rad (v*t/2/pi()/r*2*pi())
     $ende = $start_drehung + $rad_S_teilung * $sichtbar - 0.0001;
     $start_drehung = $start_drehung - $rad_S_teilung * $sichtbar;
	  # alle Schaufeln einzeln
#     $ende = 2*pi()+$start_drehung-0.0001;
     
     for ($drehung = $start_drehung; $drehung <= $ende; $drehung = $drehung + $rad_S_teilung)
     {                  # aufgrund eines rundungsfehlers wird nach der letzten Schaufel noch nicht 2pi() erreicht; daher -0.0001
     # die Schaufel # hier großer radius
     $r_km = $r_Rad - $radiale_Tiefe + $r_innen_Schaufel;
     $x_gm   = -sin($rad_Schaufel) * $r_Schaufel;
     $y_gm   = $r_Rad - cos($rad_Schaufel) * $r_Schaufel;
     $bogen  = $rad_Schaufel;
     print MTV "\n\n\%linetype = 1 linewidth = 1 linecolor = 4\n";
     do {  $xz = sin($bogen)* $r_Schaufel + $x_gm;
	   $yz = cos($bogen)* $r_Schaufel + $y_gm;

	   $x  = sqrt($xz**2 + $yz**2)* sin(atan($xz/$yz)+$drehung);
        	  # $y  = sqrt(($xz**2 + $yz**2) - $x**2); # sobald y negativ sein soll wird es wieder positiv
	   $y  = sqrt($yz**2 + $xz**2)* cos(atan($xz/$yz)+$drehung);
	   $x_km = sin($bogen)* ($r_Schaufel - $r_innen_Schaufel) + $x_gm;
	   $y_km = cos($bogen)* ($r_Schaufel - $r_innen_Schaufel) + $y_gm;
	   print MTV "$x  $y\n";
	   $bogen = $bogen + $rad_schritt/10; # 10 fache Genauigkeit sonst überschneidet die Schaufel den Boden
	}
     until $x_km**2 + $y_km**2 <= $r_km **2;

             # die Schaufel # hier kleiner Radius
	   $r_Sbq_ist = $x**2 + $y**2;  # kontrollvariable wird den letzten r^2 speichern
     for (;;)    #  for ever oder bis last
	 { $xz = sin($bogen)* $r_innen_Schaufel + $x_km;
	   $yz = cos($bogen)* $r_innen_Schaufel + $y_km;
	   $x  = sqrt($xz**2 + $yz**2)* sin(atan($xz/$yz)+$drehung);
        	    #$y  = sqrt(($xz**2 + $yz**2) - $x**2);
	   $y  = sqrt($xz**2 + $yz**2)* cos(atan($xz/$yz)+$drehung);
	   $r_Sbq_war = $r_Sbq_ist;
	   $r_Sbq_ist = $x**2 + $y**2;
	   if ( $r_Sbq_ist > $r_Sbq_war)
	    { last;
	    }
	   print MTV "$x  $y\n";
	   $bogen = $bogen + $rad_schritt;
	 }

	#print "das ist der durchgang $i\n";
	#$i++;
     }

     # Jetzt kommt die Rinne
     print MTV "\n\n%linetype = 1 linewidth = 1 linecolor = 6\n";
      $x   = $x_Rinne;
      $y   = sqrt(($s_Rinne+$r_Rad)**2 - $x**2);
      print MTV "$x  $y\n";
      print MTV (($x - cos($rad_Rinne)*0.1) . "  " . ($y + sin($rad_Rinne)*0.1). "\n"); # Einfach 10 cm lang das Ding

     # Jetzt kommt das Wasser
     print MTV "\n# das Wasser:\n";
     print MTV "\n\n%linetype = 1 linewidth = 1 linecolor = 5 linelabel = \"Wasser\"\n";
       $xz = $x;
       $yz = $y;                 # startwerte sichern

       for ($t  = $t_W ; $t < $t_W+$dt_Wasser ; $t = $t + $dt_Wasser/40)
       { $x  = $v_Wasser * $t * cos($rad_Rinne) + $xz;
	 $y  = $yz - $v_Wasser * $t * sin($rad_Rinne) - $g * $t**2 / 2 ;
	 
	 print MTV "$x  $y\n";
       }
     print MTV "\n\n%linetype = 1 linewidth = 1 linecolor = 5\n";
       $yz = $yz + cos($rad_Rinne) * $h_Wasser;         # startwert für oberkante Wasser
       $xz = $xz + sin($rad_Rinne) * $h_Wasser;

       for ($t  = $t_W ; $t < $t_W+$dt_Wasser ; $t = $t + $dt_Wasser/10)
       { 
	 $x  = $v_Wasser * $t * cos($rad_Rinne) + $xz;
	 $y  = $yz - $v_Wasser * $t * sin($rad_Rinne) - $g * $t**2 / 2 ;
	 print MTV "$x  $y\n";
       }

}
print MTV "\n\n\$END\n";

close MTV;

system  "/usr/X11R6/bin/plotmtv  -geom 660x480 -bg '#181880' -fg '#ffffff'   $mtv >/dev/null ";
    # -geom 1000x800

system "/bin/rm $mtv";