#!/usr/bin/perl

# alpha_inversion_Faran.pl
# (c) 2003-2004 Andreas Richter. All rights reserved.
# This program is free software; you can redistribute it and/or
# modify it under the same terms as Perl itself.

# alpha_inversion_Faran.pl liefert Partikelgrößenverteilungen auf der
# Basis des Modells nach Faran durch Inversion der Koeffizienten-Matrix
#
# Die Datei benötigt als Eingabe den Namen einer Textdatei, welche
# die gemessenen Extinktionsfunktionen enthält. Das Skript extrahiert
# die Daten und speichert die Mittelwerte in einem Array
#

# Schritte in der Ausführung des Skripts:
#    - Einlesen der Kommandozeilen-Parameter
#    - Überprüfen der Kommandozeilen-Parameter
#    - Einlesen von ka und Extinktionsmesswerten
#    -

# Kommandozeilenparameter
# -k <Dateiname>             : Dateiname mit den Extinktionskoeffizienten
# -h <Dateiname>             : Name der Datei, die die Partikelgrößenver-
#                              teilung enthält
# -P <Dateiname>             : Name der Datei, die die Partikelstoff-
#                              eigenschaften enthält
# -F <Dateiname>             : Name der Datei, die die Fluidstoff-
#                              eigenschaften enthält

use strict;
use warnings;

use Getopt::Std; # Komandozeilen-Optionen interpretieren

# Transformation Np/m/MHz -> dB/cm/MHz
my $dB = 20/log(10) / 100;

# Variablen

my @ka_m;          # ka-Werte (dimensionslose Wellenzahl) - gemessen
my @alpha_m;       # alpha-Werte (Dämpfungsspektrum) - gemessen
my @alpha_m_0;       # alpha-Werte (Dämpfungsspektrum) - gemessen
my $kext_file="";  # Name der Datei, die die gemessenen Extinktionskoeffizienten
                   # enthält
my $particle_file="";
my $fluid_file="";
my $result_file=""; # Name der Datei, die Ausgabe enthält
my $from_stdin=0;  # falls 1: Hole alle Werte von STDIN, ohne weitere Rückfrage
my @xo_f;          # Partikelgrößenverteilung (PGV): Klassenobergrenzen in m
my @Q3_f;          # PGV: Summenfunktion Q3 in Werten von 0 - 1
my @dQ3_f;         # PGV: dQ3 = q3*dx
my @xm_f;          # PGV: Klassenmitten in m
my @q3_f;          # PGV: Dichtefunktion q3 in 1/µm
my $x_50_ini=15;   # PGV: Medianwert in µm
my $phi_ini=0.01;  # PGV: Volumenanteil Feststoff
my $index=0;

# Frequenzen
my @f_f; #= pdl [3.0, 3.69, 4.53, 5.57, 6.84, 8.4, 10.32, 12.69, 15.59, 19.15, 23.53, 28.91, 35.53, 43.65, 53.64, 65.9, 80.9, 99.0];

my $c_F;           # Schallgeschwindigkeit im Fluid in m/s
my $c_PL;          # Schallgeschwindigkeit (longitudinal) im Partikel in m/s
my $c_PT;          # Schallgeschwindigkeit (transversal) im Partikel in m/s

my $rho_F;          # Dichte des Fluids
my $K_F;            # Kompressionsmodul Fluid
my $eta_F;          # Viskosität Fluid

my $rho_P;          # Dichte des Partikels
my $K_P;            # Kompressionsmodul Partikel

my $s=0;           # Startwerte der zu optimierenden Parameter.
my $t=0;
my $u=0;

select(STDOUT); $| = 1; # make STDOUT unbuffered


# Kommandozeilenparameter einlesen
my %options;

my @desired_opts= qw/k P F r s/;

# Parse options & store them in %options
getopt('kPFrs', \%options);

# store them in more meaningful variables
$kext_file = $options{'k'};
$particle_file = $options{'P'};
$fluid_file = $options{'F'};
$result_file = $options{'r'};
$from_stdin = $options{'s'};

# Kommandozeilenparameter überprüfen.
# Wir nehmen die Eingabe von STDIN, wenn Option "-s" angegeben wurde
# ansonsten: andere Kommandozeilenparameter auslesen und entsprechende
# Dateien importieren.
if ($result_file) {
  open(STDOUT, ">$result_file")
}
unless ($from_stdin) {
  if (!$kext_file) {
    die("Datei mit den Extinktionskoeffizienten nicht spezifiziert -
    Option -k verwenden\n");
  }

  if (!$particle_file) {
    die("Datei mit den PartikelEigenschaften nicht angegeben - Option -P verwenden.");
  }

  if (!$fluid_file) {
    die("Datei mit den Fluid-Eigenschaften nicht angegeben - Option -F verwenden.");
  }

  printf("Dämpfungsmesswerte in Datei $kext_file\n");


  # Fluid properties input
  my $linebuf;
  open(FILE,$fluid_file) or die("Can't open fluid properties file $fluid_file\r\n");
  while ( $linebuf=<FILE> ) {
    chop($linebuf);
    my @data_line=split(' ',$linebuf);
    if ($data_line[0] eq 'c') { $c_F=$data_line[1];}
    if ($data_line[0] eq 'rho') { $rho_F=$data_line[1];}
    if ($data_line[0] eq 'eta') { $eta_F=$data_line[1];}
    if ($data_line[0] eq 'K') { $K_F=$data_line[1];}
  }
  close(FILE);
  #printf("Fluid: K = $K_f, rho = $rho_f, eta = $eta_f\n");

  # Particle properties input
  open(FILE,$particle_file) or die("Can't open particle properties file $particle_file\r\n");
  while ( $linebuf=<FILE> ) {
    chop($linebuf);
    my @data_line=split(' ',$linebuf);
    if ($data_line[0] eq 'c') {$c_PL=$data_line[1];}
    if ($data_line[0] eq 'c_T') {$c_PT=$data_line[1];}
    if ($data_line[0] eq 'rho') { $rho_P=$data_line[1];}
    if ($data_line[0] eq 'K') { $K_P=$data_line[1];}
  }
  close(FILE);

  # Öffnen der Datei, die die Extinktionskoeffizienten enthält.
  open(FILE,$kext_file) or die("Kann die Dämpfungsmesswerte-Datei $kext_file nicht öffnen!\n");

  # Header lesen
  # = Zeilen einlese, bis die Zeile den gewünschten Text ("alpha_avg") enthält.
  my @fields;
  do {
    $linebuf=<FILE>;                                # get one line of input
    chomp($linebuf);
    # Die Schallgeschwindigkeit extrahieren
    if (index($linebuf,"# c=")!=-1) {
      @fields=split("\t",$linebuf);
      #$c_F=$fields[1];
    }
  } until (index($linebuf,"alpha_avg")!=-1);

  # bestimme die Position der f- und ka-Spalte und der Spalte von alpha_avg;
  @fields = split("\t", $linebuf);
  my $i=0;
  while ( ($i<=$#fields) && ($fields[$i] ne "# f [MHz]") ) {$i++;}
  my $f_pos=$i;
  #printf("f: $i\n");
  while ( ($i<=$#fields) && ($fields[$i] ne "ka") ) {$i++;}
  my $ka_pos=$i;
  #printf("ka: $i\n");
  while ( ($i<=$#fields) && ($fields[$i] ne "alpha_0") ) {$i++;}
  my $alpha_0_pos=$i;
  while ( ($i<=$#fields) && ($fields[$i] ne "alpha_avg") ) {$i++;}
  my $alpha_pos=$i;
  #printf("alpha: $i\n");

  # Datenzeilen einlesen und Daten (f, ka, alpha) extrahieren.
  $i=0;
  do {
    $linebuf=<FILE>;
    chomp($linebuf);
    @fields = split("\t", $linebuf);
    #printf("@fields\n");
    unless ($fields[$alpha_pos] eq "-") {
      $alpha_m_0[$i]=$fields[$alpha_0_pos];
      $alpha_m[$i]=$fields[$alpha_pos];
      $ka_m[$i]=$fields[$ka_pos];
      $f_f[$i] = $fields[$f_pos];
      $i++;
    }
  } until ( eof(FILE));

  close(FILE);
}
else {
  # Option -s wurde angegeben - STDIN auslesen /parsen.
  my $linebuf;
  my $i=0;
  # 1. Die Dämpfungswerte lesen.
  while (defined($linebuf=<STDIN>)) {
    chomp($linebuf);
    if (($linebuf=~"Begin of attenuation data")) {
    # 1. Die Dämpfungswerte lesen.
      while (defined($linebuf=<STDIN>) && !($linebuf=~"End of attenuation data")) {
      # kommt in Spalten, tab-Separiert
      # Frequenz \t Dämpfung des reinen Fluids \t Dämpfung der Suspension
        chomp($linebuf);
        my @fields=split("\t",$linebuf);
        unless ($fields[2] eq "-") {
          $f_f[$i]=$fields[0];
          $alpha_m_0[$i]=$fields[1];
          $alpha_m[$i]=$fields[2];
          $i++;
        }
      }
    }
    if (($linebuf=~"Begin of particle properties data")) {
    # 2. Partikeleigenschaften lesen.
      while (defined($linebuf=<STDIN>) && !($linebuf=~"End of particle properties data")) {
        # Als Zeilen, mit Erkennungswörten "c", "c_T" und "rho" beginnend
        # Zelle nach dem darauffolgenden Leerzeichen beinhaltet Wert.
        chomp($linebuf);
        my @data_line=split(' ',$linebuf);
        if ($data_line[0] eq 'c') {$c_PL=$data_line[1];}
        if ($data_line[0] eq 'c_T') {$c_PT=$data_line[1];}
        if ($data_line[0] eq 'rho') { $rho_P=$data_line[1];}
      }
    }
    if (($linebuf=~"Begin of fluid properties data")) {
    # 3. Eigenschaften der kontinuierlichen Phase lesen.
      while (defined($linebuf=<STDIN>) && !($linebuf=~"End of fluid properties data")) {
        # Als Zeilen, mit Erkennungswörten "c", und "rho" beginnend
        # Zelle nach dem darauffolgenden Leerzeichen beinhaltet Wert.
        chomp($linebuf);
        my @data_line=split(' ',$linebuf);
        if ($data_line[0] eq 'c') { $c_F=$data_line[1];}
        if ($data_line[0] eq 'rho') { $rho_F=$data_line[1];}
      }
    }
    if (($linebuf=~"Begin of initial data")) {
      while (defined($linebuf=<STDIN>) && !($linebuf=~"End of initial data")) {
        # Als Zeilen, mit Erkennungswörten "c", und "rho" beginnend
        # Zelle nach dem darauffolgenden Leerzeichen beinhaltet Wert.
        chomp($linebuf);
        my @data_line=split(' ',$linebuf);
        if ($data_line[0] eq 'x_50') { $x_50_ini=$data_line[1];}
        if ($data_line[0] eq 'phi') { $phi_ini=$data_line[1];}
      }
    }
  }
    # Hier sollte ein Test erfolgen, der überprüft, ob alle Daten korrekt eingelesen wurden.
}

# Berechne Dichte- und Geschwindigkeitsparameter
my $g = $rho_P / $rho_F;
my $h_L = $c_PL / $c_F;
my $h_T = $c_PT / $c_F;


#printf("@alpha_m\n @ka_m\n @f_f\n");

# Kurvenanpassung starten

# Pfad zu den Modellen spezifizieren
use lib qw( frontend/Module Module );

use PDL;
use PDL::Math;
use PDL::Primitive;
use PDL::Complex;

# Modul mit dem funktionellen Zusammenhang zwischen x, f und alpha
use Faran_lsg_pdl;
use PKM;
use Helos;


# Dämpfungsspektrum + Frequenzen -> PDL
my $alpha = pdl @alpha_m;
my $f = pdl @f_f;


# Partikelgrößenklassen ermitteln
my $x;
my $x50;

($x,$x50) = calculate_x($f,$alpha,$x_50_ini);

$x*=1.5;
$x->slice(0).=0.005e-6;
# Startwert für die Volumenkonzentration
my $phi=$phi_ini;
#-----------------------------------------------
# Lösungsverfahren:
#-----------------------------------------------

use PDL::Slatec;
my $dQ;

my $K=alpha_Str_Faran_noweight($f, $x, $c_F, $h_L, $h_T, $g, $phi)*$dB;


# 0. Condition number
my ($lsg,$condition_number,$z) = geco($K);

printf "Condition number $condition_number\t $z\t $lsg\n";
# 1. Einfache Inversion von $K

#$dQ = matmult(matinv($K),transpose($alpha));

# 2. Phillips - Twomey nach Riebel.

# Bilde die Glättungsmatrix
my $H = zeroes($alpha->getdim("0"), $alpha->getdim("0"));
$H->diagonal(0,1).=-1;


# in den PDLs sind die Indices vertauscht (im Vergleich zu
# herkömmlicher Matrix-Notation.
# PDL: zuerst Spalte, dann Zeile
$H->slice("0,0").=0;
$H->slice("1,1").=0;

for(my $i=2;$i<$H->getdim("0");$i++){
  my $j=$i-1;
  my $k=$j-1;
  $H->slice("$j,$i").=2;
  $H->slice("$k,$i").=-1;
}

my $gamma=0.001; # Glättungsfaktor
# $alpha muss transponiert werden, weil die Variable als Zeilenvektor
# gespeichert wird.
#$dQ = matmult(matmult(matinv(matmult(transpose($K),$K)+ $gamma* matmult(transpose($H),$H) ),transpose($K)),transpose($alpha));

# 3. P-T mit modifizierte Glättungsmatrix (laut Kräuter, Diss., S.15)
my $H_opt = matmult(matinv(matmult(transpose($K),$K)),transpose($K));
my $peaks;
# do {
#   $gamma*=10;
#   $dQ = matmult(matmult(matinv(matmult(transpose($K),$K)+ $gamma* matmult(transpose($H_opt),$H_opt) ),transpose($K)),transpose($alpha));
#   $peaks=no_of_peaks(transpose($dQ));
#   printf $peaks . "\t" . $gamma . "\n";
# } until ($peaks<=1);



# 4. P-T mit modifizierte Glättungsmatrix (laut Kräuter, Diss., S.15)
# + Anpassung des Konzentrationsparameters

# Schwellwert (Abbruchkriterium) der Konzentrationsmessung
my $epsilon=0.0001;
my $phi_prev;
my $MSQ=10000;
my $MSQ_old;
do {
#  $phi=0.011;
#  $phi_prev=0.011;
  $K=alpha_Str_Faran_noweight($f, $x, $c_F, $h_L, $h_T, $g, $phi)*$dB;
  unless ($from_stdin) {
    # calculate viscoinertial attenuations
    #richtig
    my $res_vis_noweigh =alpha_visc_noweight($f,$x/2,$phi,$K_P, $K_F, $rho_P, $rho_F, $eta_F);
    # falsch, aber an Prof. Klöden übermittelt.
    #my $res_vis_noweigh =alpha_visc_noweight($f,$x/2,$phi/100,$K_P, $K_F, $rho_P, $rho_F, $eta_F);
    $K+=$res_vis_noweigh;
  }
  printf $K . "\n";
  # printf "K ok \n";
  $H_opt = matmult(matinv(matmult(transpose($K),$K)),transpose($K));
  # printf "H_opt ok \n";
  $gamma=0.001;
  # Adapt gamma
  do {
    $gamma*=2;
    $dQ = matmult(matmult(matinv(matmult(transpose($K),$K)+ $gamma* matmult(transpose($H_opt),$H_opt) ),transpose($K)),transpose($alpha));
    $peaks=no_of_peaks(transpose($dQ));
  } until ($peaks==2)|| ($gamma >= 1e270);
  #$dQ = matmult(matmult(matinv(matmult(transpose($K),$K)+ $gamma* matmult(transpose($H_opt),$H_opt) ),transpose($K)),transpose($alpha));
  #printf sumover(transpose($dQ)) . "\n";
  # Setze alle negativen Werte auf 0
  (my $dQ_0=$dQ->where($dQ<0)) .=0;
  #printf $dQ_0 . "\n";
  # normalisiere verbliebende Werte
  $dQ=$dQ/sumover(transpose($dQ));
  # printf "dQ ok \n";
  my $alpha_res=alpha_Str_Faran($f, $x, transpose($dQ), $c_F, $h_L, $h_T, $g, $phi)*$dB;
  unless ($from_stdin) {
    # calculate viscoinertial attenuations
    # richtig
    my $res_vis_noweigh =alpha_visc($f,$x/2,transpose($dQ),$phi,$K_P, $K_F, $rho_P, $rho_F, $eta_F);
    #printf "res_vis\n";
    # falsch, aber an Prof. Klöden übermittelt
    #my $res_vis_noweigh =alpha_visc($f,$x/2,transpose($dQ),$phi/100,$K_P, $K_F, $rho_P, $rho_F, $eta_F);
    $alpha_res+=$res_vis_noweigh;
  }
  # printf "alpha_Str ok \n";
  $phi_prev=$phi;
  # printf "phi_prev ok \n";
  $phi=$phi_prev * ($alpha->max/$alpha_res->max);
  # printf "phi ok \n";
  # Abbruchkriterium:
  $MSQ_old=$MSQ;
  $MSQ=sumover(($alpha_res-$alpha)**2);
  printf "phi " . $phi . " MSQ ". $MSQ ."\n";
#  my $a=graph_dQ($x,$dQ);
}
until ((abs($phi_prev-$phi)<$epsilon)|| ($MSQ_old<=$MSQ));

$K=alpha_Str_Faran_noweight($f, $x, $c_F, $h_L, $h_T, $g, $phi)*$dB;
unless ($from_stdin) {
  # calculate viscoinertial attenuations
  # richtig
  my $res_vis_noweigh =alpha_visc_noweight($f,$x/2,$phi,$K_P, $K_F, $rho_P, $rho_F, $eta_F);
  # falsch, aber an Prof. Klöden übermittelt
  #my $res_vis_noweigh =alpha_visc_noweight($f,$x/2,$phi/100,$K_P, $K_F, $rho_P, $rho_F, $eta_F);
  $K+=$res_vis_noweigh;
}


printf "c_F " . $c_F . "\n";
printf "h_L " . $h_L . "\n";
printf "h_T " . $h_T . "\n";
printf "g " . $g . "\n";
printf "phi " . $phi . "\n";
printf "epsilon " . $epsilon . "\n";
printf "gamma " . $gamma . "\n";
printf "Kern-Matrix:\n";
printf $K . "\n";

printf "Messwerte-Vektor:\n";
printf $alpha . "\n";

# Calculate dq*
my $dx = calculate_dx($x);
my $q = transpose($dQ)->clump(2)/$dx;
my $xm = xm_from_xo($x);

#-----------------------------------------------
# Ausgabe:
#-----------------------------------------------

# # Differenz (Inversionsergebnis - Messung)
# # 1. Rückrechnung durch Multiplikation von $dQ mit der Matrix $K
# #    Abziehen von $alpha liefert den Fehler der Methode
# my $res=inner($K,transpose($dQ))-$alpha;
#
# # 2. Rückrechung durch direkte Bildung des Dämpfungsspektrums
# # die Verknüpfung mit $dQ geschieht in alpha_Str_Faran
# my $alpha_res=alpha_Str_Faran($f, $x, transpose($dQ), $c_F, $h_L, $h_T, $g, $phi)*$dB;
# my $res_test=$alpha_res-$alpha;
#
# # Ergebnisse sollten identisch sein.
# printf("res $res\n");
# printf("res_test $res_test\n");

#printf("H $H\n");
#printf("H_opt $H_opt\n");

printf "x " . $x ."\n";
printf "dQ " . transpose($dQ)->clump(2) . "\n";

# activate this for scripted calculation and graph creation
print STDOUT ("xo\tdQ_fit\txm\tq_fit\n");
for (my $i=0; $i<$x->getdim("0");$i++) {
  my $Q_out = transpose($dQ)->clump(2)->slice("$i")->sclr;
  my $xo = $x->slice("$i")->sclr;
  my $xm_out = $xm->slice("$i")->sclr;
  my $q_out = $q->slice("$i")->sclr; # Unit: 1/m
  print STDOUT ("$xo\t$Q_out\t$xm_out\t$q_out\n");
}


#printf cumusumover(transpose($dQ)) . "\n";


# Graphische Darstellung
# pgplot
#graph_dQ($x,$dQ);
# gnuplot->eps
#plot_PSD_graph($x,cumusumover(transpose($dQ)->clump(2)),$xm,$q);


# Darstellung des resultierenden Dämpfungsspektrums
# my $alpha_res=alpha_Str_Faran($f, $x, transpose($dQ), $c_F, $h_L, $h_T, $g, $phi)*$dB;
# dev('/XSERVE');
# line($f,$alpha, {LINESTYLE=>1});
# hold;
# points($f,inner($K,transpose($dQ)), {SYMBOL=>1});
# points($f,$alpha_res, {SYMBOL=>2});


exit(0);


# ----------------------------------------
# Subs
# ----------------------------------------

# Calculate the number of peaks a curve has
sub no_of_peaks {
  # vector of values which contain maxima and minima
  my ($y)=@_;
  my $no_values=$y->getdim(0);
  # minimum of vector
  my $min_y=$y->min();
  # work with a copy of the original vector
  # (because values get deleted in the process)
  my $test_vector = $y->copy;
  # initialize a mask that shows which values have already
  # been identified as "current maximum"
  my $mask=zeroes($test_vector->getdim(0));
  # initialize the number of peaks
  my $no_of_peaks=0;
  # test all except one value (the minimum)
  for(my $i=0;$i<$no_values-1;$i++) {
    # get index of current maximum
    my $max_index=maximum_ind($test_vector)->sclr;
    # unless the value is not at the border of the vector
    unless (($max_index<=0) || ($max_index>=$no_values-1)) {
      # unless the neighboring values have already been found as peaks or parts thereof
      unless ($mask->slice($max_index+1) || $mask->slice($max_index-1)) {
        # it is a peak. increase number of peaks
        $no_of_peaks++;
      }
    # if value is at the border of the vector,
    } else {
      # test if it is large enough to count as a separate peak
      # (adjust the threshold to your needs)
      if ($test_vector->slice($max_index)->sclr>0.0001)  {
        # increase number of peaks
        $no_of_peaks++;
      }
    }
    # set value to be equal to the minimum, just so that we don't find it again
    $test_vector->slice($max_index).=$min_y;
    # mark the value as "alread found" (=1 in the mask)
    $mask->slice($max_index).=1;
  }
  # return the number of peaks
  return $no_of_peaks;
}


# # Create a PGPLOT graph that shows the particle size distributions
# sub graph_dQ {
#   use PDL::Graphics::PGPLOT;
#   use PGPLOT;
# 
#   my ($x,$dQ)=@_;
#   # Maxima und Minima der Darstellung
#   my $ymax = $dQ->max;
#   my $ymin = $dQ->min;
#   if ($ymin<0) {$ymin=$ymin*1.1}
#     else {$ymin=$ymin*0.9};
#   if ($ymax<0) {$ymax=$ymax*0.9}
#     else {$ymax=$ymax*1.1};
#   my $xmax =$x->max *1e6;
#   my $xmin =$x->min *1e6;
# 
#   dev('/XSERVE');
#   #dev("Phillips_Twomey_mod_gamma$gamma.eps\/PS");
# 
#   # Darstellung der PGV-Messung (Vergleich mit HELOS)
#   my $bord_opt = { TYPE => 'RELATIVE', VALUE => 0.1 };
#   env($xmin,$xmax,$ymin, $ymax, {YRANGE=>[$ymin,$ymax], BORDER=>$bord_opt});
#   line($x*1e6,$dQ, {SYMBOL=>2,LINESTYLE=>1});
#   hold;
#   points($x*1e6,$dQ, {SYMBOL=>2});
# 
#   label_axes('x [µm]', 'dQ_3', '');
#   #line($x*1e6,$res, {LINESTYLE=>2});
#   legend(['US-Inversion'], $xmax*0.6, $ymax*0.8, { symbol =>[2]});
#   close_window();
# }

# Führt eine Berechnung der Stützstellen der PGV durch.

# Der Mittelpunkt des (logarithmierten) Vektors der Stütztellen
# wird anhand einer empirischen Funktion bestimmt.
# Diese Funktion ordnet dem Maximum des Dämpfungsspektrums eine
# Partikelgröße zu, und wird als Mittelpunkt des Vektors gesetzt.
# Ziel: Stützstellen liegen im Bereich der zu berechnenden PGV

sub calculate_x {
  my ($f, $alpha,$x_50_ini) = @_;
  # Parameter der empririschen Funktion x = f(f_max)
  # Ermittelt für Glaspartikel in Ethanol, Wasser und NaCl-Lösung
  my $n = -0.98475052;
  my $K = 3.0542539;
  # Bestimme Vektor-Index des Maximums des Dämpfungsspektrums
  my $fmax_index = maximum_ind($alpha);
  # Bestimme zugehörige Frequenz
  my $fmax = $f->slice("$fmax_index");
  #print STDOUT  "fmax " . $fmax . "\n";
  # Bestimme zugehörige Frequenz
  my $x50;
  if ($x_50_ini>0) {
    $x50=$x_50_ini;
  } else {
    $x50=10**$K * ($fmax**$n);
  }
  #print STDOUT  "x50 " . $x50 . "\n";
  # Bestimme Parameter logarithmischen Reihe
  #my $dx_ln = log(1/0.90);
  my $dx_ln = log(1/0.85);
  #my $dx_ln = log(1/0.80);
  #my $dx_ln = log(1/0.78);
  #my $dx_ln = log(1/0.75);
  # Länge des Partikelgrößenvektors definiert sich aus Anzahl der Frequenzen
  my $x=sequence($f->getdim("0"));
  #print STDOUT  "x " . $x . "\n";
  # Anzahl der Frequenzen, geteilt durch 2, ergibt den Mittelwert (Index)
  my $a=floor(($f->getdim("0")) / 2);
  #$a=rint(8.5);
  #print STDOUT  "a " . $a . "\n";
  $x = $x50 * exp(($x - $x->slice("$a"))*$dx_ln);  
  #$x=$x->rotate(1);
  # Rückgabe. x = Klassenmittelwerte. Rückgabewerte sind in m
  print STDOUT  "x " . $x . "\n";
  return ($x/1e6, $x50/1e6);
}


sub plot_PSD_graph {
  my ($xo,$Q,$xm,$q)=@_;
  my $gnuplot = '/usr/bin/gnuplot';
  my $boundary1=min($xm)*0.8*1e6;
  my $boundary2=max($xo)*2*1e6;

  # grün -> dunkelrot
  # blau -> dunkelgrün
  # violett -> dunkelblau
  # türkis -> dunkelgelb
  # gelb -> dunkeltürkis
  # open ( GNUPLOT, "|$gnuplot | sed 's/0 1 0/0.752941 0 0/;
  #                                   s/0 0 1/0 0.752941 0/;
  #                                   s/1 0 1/0 0 0.752941/;
  #                                   s/0 1 1/0.752941 0.752941 0/;
  #                                   s/1 1 0/0 0.752941 0.752941/'");

  #                                   R G B R  G  B
  open ( GNUPLOT, "|$gnuplot | sed \""
                                  .'s/PL \[\] 1 0 0/PL \[2 dl 4 dl\]0.752941 0 0/;' # 1 rot -> dunkelrot
                                  .'s/PL \[4 dl 2 dl\] 0 1 0/PL \[2 dl 4 dl\] 0 0.752941 0/g;' # 2 grün -> dunkelgrün
                                  .'s/PL \[2 dl 3 dl\] 0 0 1/PL \[2 dl 4 dl\] 0 0 0.752941/;' # 3 blau -> dunkelblau
                                  .'s/PL \[1 dl 1.5 dl\] 1 0 1/PL \[\] 0.752941 0 0/;' # 4 magenta ->dunkelrot
                                  .'s/PL \[5 dl 2 dl 1 dl 2 dl\] 0 1 1/PL \[\] 0 0.752941 0/;' # 5 cyan -> dunkelgrün
                                  .'s/PL \[4 dl 3 dl 1 dl 3 dl\] 1 1 0/PL \[\] 0 0 0.752941/;' # 6 gelb -> dunkelblau
                                  .'s/PL \[2 dl 2 dl 2 dl 4 dl\] 0 0 0/PL \[\] 0 0 0/;' # 7 schwarz -> dunkelgrün
                                #  .'s/PL \[2 dl 2 dl 2 dl 2 dl 2 dl 4 dl\] 1 0.3 0/PL \[\] 0 0 0.6/;' # 8 orange -> dunkelblau
                                #  .'s/PL \[2 dl 2 dl 2 dl 2 dl 2 dl 2 dl 2 dl 4 dl\] 0.5 0.5 0.5/PL \[\] 0.6 0.6 0/;'# 9 grau ->  dunkelgelb
                                  ."\" > Glas_H2O_40_50_1Vol_PT_gamma100.eps" );

  #open (GNUPLOT,">test.txt");
  print GNUPLOT "\n";
  print GNUPLOT << "gnuplot_Commands";
  # color plor
  set terminal postscript eps enhance color solid "Helvetica" 23
  # 23
  set size 0.857, 0.680
  set missing "-"
  set ytics nomirror
  set xtics nomirror
  set y2tics
  set logscale x
  set xrange [$boundary1:$boundary2]
  set yrange [0:1]
  set y2range [0:]
  set xlabel "x_{o} [{\/Symbol m}m]"
  set ylabel "Q_{3}"
  set y2label "q_{3}^{*}"
  #rot
  set linestyle 1 linetype 1 lw 4 pointtype 7 pointsize 0.75
  #grün
  set linestyle 2 linetype 2 lw 4 pointtype 6 pointsize 0.75
  #blau
  set linestyle 3 linetype 3 lw 4 pointtype 9 pointsize 0.75
  #magenta
  set linestyle 4 linetype 4 lw 4 pointtype 8 pointsize 0.75
  #cyan
  set linestyle 5 linetype 5 lw 4 pointtype 5 pointsize 0.75
  #gelb
  set linestyle 6 linetype 6 lw 4 pointtype 7 pointsize 0.75
  #schwarz
  set linestyle 7 linetype 7 lw 4 pointtype 6 pointsize 0.75
  #orange
  set linestyle 8 linetype 8 lw 4 pointtype 9 pointsize 0.75
  #grau
  set linestyle 9 linetype 9 lw 4 pointtype 8 pointsize 0.75
  #rot
  set linestyle 10 linetype 10 lw 4 pointtype 5 pointsize 0.75
  set key Left left reverse
  set key spacing 1.3
gnuplot_Commands
  my $q_c_star=zeroes($q->getdim(0));
  $q_c_star.=$xm/log10(log(10))*$q;
  my $gnuplotcommand_string ="plot ";
  my $gnuplotdata;
  # Q3 on first axis
  $gnuplotcommand_string .= "\'-\' using 1:2 notitle with lines linestyle ";
  $gnuplotcommand_string .= 1;
  $gnuplotcommand_string .= ", ";
  # q3 - on 2nd axis
  $gnuplotcommand_string .= "\'-\' using 1:2 axes x1y2 notitle with lines linestyle ";
  $gnuplotcommand_string .= 3;
  $gnuplotcommand_string .= "\n";
    # Q3 Data
  $gnuplotdata .= (2*$xo->slice("0")-$xo->slice("1"))->sclr*1e6 . "\t0\n";
  for(my $i=0; $i<$Q->getdim(0);$i++) {
    $gnuplotdata .= $xo->slice("$i")->sclr*1e6 ."\t". $Q->slice("$i")->sclr ."\n";
  }
  $gnuplotdata .= "e\n";
  # q3 Data
  for(my $i=0; $i<$q->getdim(0);$i++) {
    $gnuplotdata .= $xm->slice("$i")->sclr*1e6 ."\t". $q_c_star->slice("$i")->sclr ."\n";
  }
  $gnuplotdata .= "e\n";
  print GNUPLOT $gnuplotcommand_string;
  #print STDOUT $gnuplotcommand_string;
  print GNUPLOT $gnuplotdata;
  #print STDOUT $gnuplotdata;
  close (GNUPLOT);
}
