#!/usr/bin/perl

# kext_schaetz_LNVT_Faran.pl
# (c) 2004-2005 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_schaetz_LNVT_Faran_PGV_calc.pl liefert Partikelgrößenverteilungen
# auf der Basis des Modells nach Faran durch Schätzung der Parameter
# einer logarithmischen Normalverteilung
#
# 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
# -P <Dateiname>             : Name der Datei, die die Partikelstoff-
#                              eigenschaften enthält
# -F <Dateiname>             : Name der Datei, die die Fluidstoff-
#                              eigenschaften enthält
# -r <Dateiname>             : Name der Datei, in der die Ergebnisse gespei-
#                              chert werden
# -s                         : wenn angegeben, werden alle anderen Optionen
#                              ignoriert, und die Daten von STDIN gelesen
# -S mu,sigmaL-ln,phi        : Startwerte für die Anpassparameter


use strict;
use warnings;

use Getopt::Std; # Komandozeilen-Optionen interpretieren
# Variablen

my @ka_m;          # ka-Werte (dimensionslose Wellenzahl) - gemessen
my @kext_m;        # kext-Werte (Extinktionskoeffizienten) - gemessen
my @alpha_m;       # Dämpfungswerte - gemessen
my @alpha_m_0;      # Dämpfungswerte - Nullmessung
my $kext_file="";  # Name der Datei, die die gemessenen Extinktionskoeffizienten
                   # enthält
my $alpha0_file="";  # Name der Datei, die die gemessene Dämpfung der Leermessung
                   # enthält
my $result_file=""; # Name der Datei, die Ausgabe enthält
my $particle_file=""; # Enthält Stoffwerte des Partikels
my $fluid_file=""; # Enthält Stoffwerte des Fluids
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 $index=0;

my @start_values;  # Startwerte (mu,sigma_ln,phi) von der Kommandozeile

# 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;           # Dichte des Fluids
my $rho_P;         # Dichte des Partikels

my $K_P;
my $K_F;
my $eta_F;

#my $x_50_ini=90/1e6;# Startwerte der zu optimierenden Parameter.
#my $sigma_ln_ini=0.4;
#my $phi_ini=0.02;
my $x_50_ini=84/1e6;# Startwerte der zu optimierenden Parameter.
my $sigma_ln_ini=0.42;
my $phi_ini=0.04;

my $alpha_abs_PL1;
my $alpha_abs_PL2;
my $alpha_abs_PT1;

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

# Kommandozeilenparameter einlesen
my %options;

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

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

# Parse options & store them in %options
getopt('kPFrsS', \%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 Dämpfungswerten nicht spezifiziert - Option -k verwenden\n");
  }

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

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

#  if ($result_file) {
#    open(STDOUT, ">$result_file")
#  }
  print STDOUT ("Extinktionsfunktion in Datei $kext_file\n");

  if ($options{'S'}) {
    @start_values=split(',',$options{'S'});
    unless ($#start_values==2) {
      die ("Zahl der übergebenen Startwerte stimmt nicht (mu,sigma_ln, phi)\n");
    }
    $x_50_ini=$start_values[0];
    $sigma_ln_ini=$start_values[1];
    $phi_ini=$start_values[2];
  }
  # 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=$data_line[1];}
  if ($data_line[0] eq 'K') { $K_F=$data_line[1];}
  if ($data_line[0] eq 'eta') { $eta_F=$data_line[1];}
  }
  close(FILE);
  #print STDOUT ("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/Dämpfungsspektren enthält.
  open(FILE,$kext_file) or die("Kann die Messwert-Datei $kext_file nicht öffnen!\n");

  # Header lesen
  # = Zeilen einlese, bis die Zeile den gewünschten Text ("Kext_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 ka-Spalte und der Spalte von KExt_avg;
  @fields = split("\t", $linebuf);
  my $i=0;
  while ( ($i<=$#fields) && ($fields[$i] ne "# f [MHz]") ) {$i++;}
  my $f_pos=$i;
  #print STDOUT ("f: $i\n");
  while ( ($i<=$#fields) && ($fields[$i] ne "ka") ) {$i++;}
  my $ka_pos=$i;
  #print STDOUT ("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;

  while ( ($i<=$#fields) && ($fields[$i] ne "KExt_avg") ) {$i++;}
  my $KExt_pos=$i;
  #print STDOUT ("KExt: $i\n");

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

  close(FILE);

  #print STDOUT ("@kext_m\n @ka_m\n @f_f\n");
}
else {
  # Option -s wurde angegeben - STDIN auslesen /parsen.
  my $linebuf;
  my $i=0;
  # 1. Die Dämpfungswerte lesen.
  while (defined($linebuf=<STDIN>)) {
    chomp($linebuf);
    #printf $linebuf . "\n";
    if ($linebuf=~"Begin of attenuation data") {
    # 1. Die Dämpfungswerte lesen.
      while (defined(chomp($linebuf=<STDIN>)) && !($linebuf=~"End of attenuation data")) {
      #printf $linebuf . "\n";
      # kommt in Spalten, tab-Separiert
      # Frequenz \t Dämpfung des reinen Fluids \t Dämpfung der Suspension
        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(chomp($linebuf=<STDIN>)) && !($linebuf=~"End of particle properties data")) {
      #printf $linebuf . "\n";
        # Als Zeilen, mit Erkennungswörten "c", "c_T", "rho", "K_P" 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 ($data_line[0] eq 'K') { $K_P=$data_line[1];}
      }
    }
    if (($linebuf=~"Begin of fluid properties data")) {
    # 3. Eigenschaften der kontinuierlichen Phase lesen.
      while (defined(chomp($linebuf=<STDIN>)) && !($linebuf=~"End of fluid properties data")) {
      #printf $linebuf . "\n";
        # Als Zeilen, mit Erkennungswörten "c","rho", "K_F", "eta_F" 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=$data_line[1];}
        if ($data_line[0] eq 'K') { $K_F=$data_line[1];}
        if ($data_line[0] eq 'eta') { $eta_F=$data_line[1];}
      }
    }
    if (($linebuf=~"Begin of initial data")) {
      while (defined(chomp($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 'sigma_ln') { $sigma_ln_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;
my $h_L = $c_PL/ $c_F;
my $h_T = $c_PT/ $c_F;

# Pfad zu den Modellen spezifizieren
# Wenn von Verzeichnis weiter oben getstartet:
#use lib qw( frontend/Module Module );
# Wenn vom Frontend aus gestarted
use lib qw( Module );
# Funktioniert immer auf MVT147
#use lib qw( /home/andreas/Ultraschallspektroskopie/kext_fit/Module/ );
use Helos;

use PDL;
use PDL::Core;
use PDL::Math;
use PDL::Primitive;
use PDL::Complex;
use PDL::Fit::LM; # Levenberg-Marquardt-Fit
#use PDL::Graphics2D;
#use PDL::Graphics::PGPLOT;

#$alpha_abs_PL1= (pdl[0.01438141, 0.021757636, 0.032791054, 0.049575757, 0.074760323, 0.11275025, 0.17018385, 0.25732512, 0.3883749, 0.58599843, 0.88471361, 1.3355341, 2.0172007, 3.0445776, 4.5976497, 6.9395106, 10.458143, 15.661278]);
#$alpha_abs_PL2=(pdl[ 1649626.1, 1829509.1, 2027058.2, 2247705.6, 2490768.8, 2760177.8, 3059338.3, 3392394.9, 3759965.8, 4167034.6, 4618823.8, 5119376.8, 5674883.7, 6289413.8, 6971265.1, 7725872.6, 8558601, 9465728.4]);
#$alpha_abs_PT1=(pdl[0.047981413, 0.07259108,0.10940242,0.16540206, 0.24942657, 0.37617427, 0.56779282, 0.85851602, 1.2957415, 1.9550797, 2.9516942, 4.4557841, 6.7300562, 10.157737, 15.339317, 23.152566, 34.891925, 52.251353]);


# Kurvenanpassung starten

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

my $f = pdl @f_f;
my $kext = (pdl @kext_m);
my $alpha = (pdl @alpha_m);
my $alpha_0 = (pdl @alpha_m_0);

$alpha_abs_PL1=zeroes($f) ;
$alpha_abs_PL2=zeroes($f);
$alpha_abs_PT1=zeroes($f);


my ($xm,$x50) = calculate_x($f,$alpha);
$xm*=4;
my $dx = differentiate($xm);

my $no_of_x=$xm->getdim("0");
my $dQ_result = zeroes($no_of_x);
#print STDOUT ("dQ_result $dQ_result\n");

#my $mu;

# Kurvenanpassung LNVT an die gemessenene d_kext - Kurve
my $params_LNVT_test= pdl [0,0,0];
#x50 [3.6789565e-05], sigma_ln 0.202979116544009, c [ 1.3270794]
#$params_LNVT_test->slice("0") .= 3.6789565e-05;
$params_LNVT_test->slice("0") .= $x_50_ini; #$x50; #
#$params_LNVT_test->slice("1") .= 0.302979116544009;
$params_LNVT_test->slice("1") .= $sigma_ln_ini;
$params_LNVT_test->slice("2") .= $phi_ini;

my $dQ_fit = dQ_lognorm_from_parameters_2($xm, $params_LNVT_test->slice("0"), $params_LNVT_test->slice("1"));

my $alpha_fit = alpha_ges_Faran($f,$alpha_abs_PL1, $alpha_abs_PL2, $alpha_abs_PT1,
     $xm, $dQ_fit, $c_F, $h_L, $h_T, $g, $params_LNVT_test->slice("2"))*$dB*($phi_ini);
 #$alpha_fit+=alpha_visc($f,
 #                          $xm/2,
 #                          $dQ_fit,
 #                          $phi_ini, $K_P, $K_F, $rho_P, $rho, $eta_F);

 $alpha_fit+= (1-$phi_ini)*$alpha_0;

unless ($from_stdin) {
  
  #$params_LNVT_test->slice("2") .= max($alpha_fit)/max($alpha);
}
#$alpha_fit = alpha_ges_Faran($f, $alpha_abs_PL1, $alpha_abs_PL2, $alpha_abs_PT1,$xm, $dQ_fit, $c_F, $h_L,
#             $h_T, $g, $params_LNVT_test->slice("2"))*$dB
#            +(1-$params_LNVT_test->slice("2"))*$alpha_0 ;
#printf "alpha_0 " . $alpha_0 . "\n";
unless ($from_stdin) {
  dev('/XSERVE');
}
# Startwerte
#points($f,$alpha, {SYMBOL=>4,YRANGE=>[0,max($alpha)*1.2]}); #
#hold;
#line($f,$alpha); #
#points($f,$alpha_fit,{SYMBOL=>3});
# Weighting of the data.
my $sig = (log($f)/ log($f->slice("8")))**10;
#my $sig = 1;
#my $parameter_LNVT_neu=$params_LNVT_test;
#print STDOUT "Starting fit process\n";
my ($alpha_test, $parameter_LNVT_neu, $covar_LNVT, $iterations_LNVT) =
  lmfit $f,$alpha, $sig, \&alpha_LNVT, $params_LNVT_test, {Maxiter => 50, Eps => 1e-1};

# Calculate Results
my $x50_fit=$parameter_LNVT_neu->slice("0");
my $sigma_ln_fit = $parameter_LNVT_neu->slice("1");
my $phi_fit = $parameter_LNVT_neu->slice("2");
$dQ_fit = dQ_lognorm_from_parameters_2($xm, $x50_fit, $sigma_ln_fit);
$alpha_fit = alpha_ges_Faran($f,  $alpha_abs_PL1, $alpha_abs_PL2, $alpha_abs_PT1,
     $xm, $dQ_fit, $c_F, $h_L, $h_T, $g, $phi_fit)*$dB
     +(1-$phi_fit)*$alpha_0;
#$alpha_fit+=alpha_visc($f,
#                           $xm/2,
#                           $dQ_fit,
#                           $phi_fit, $K_P, $K_F, $rho_P, $rho, $eta_F);
print STDOUT "x50 " .$x50_fit->sclr . ", sigma_ln " . $sigma_ln_fit->sclr . ", phi " . $phi_fit->sclr . "\n";

# Residuen:
#my $alpha_resid = $alpha-$alpha_fit;

# Graphische Darstellung:
# Auskommentieren für Windows
#release;
# Startwerte
#unless ($from_stdin) {
#  points($f,$alpha, {SYMBOL=>4,YRANGE=>[0,max($alpha)*1.2]}); #
#  hold;
#  line($f,$alpha); #
#  points($f,$alpha_fit,{SYMBOL=>3});
#}
#my $B = Bestimmtheitsmass($alpha,$alpha_fit);
#print STDOUT ("B $B\n");
my $s_model=s_2($alpha,$alpha_fit);
print STDOUT "s2 ". $s_model ."\n";

#Ausgabe der Ergebnisse
# unkommentieren, falls Analyse mit den Skripten notwendig ist.
#print STDOUT ("xm\tdQ_fit\n");
# Nächste Zeile für Analyse mit frontend erforderlich
print STDOUT ("xo\tdQ_fit\txm\tq_fit\n");
for (my $i=0; $i<$xm->getdim("0");$i++) {
  my $b = $dQ_fit->slice("$i")->sclr;
  my $xm = $xm->slice("$i")->sclr;
  print STDOUT ("$xm\t$b\n");
}

exit(0);

sub s_2 {
  my ($meas, $model) = @_;
  # check dimensions
  unless ($meas->getdim(0)==$model->getdim(0)) {
    return -1;
  }
  # calculate
  my $s_model2=sumover(($model-$meas)**2)/($meas->getdim("0")-1);
  return $s_model2;
}

sub Bestimmtheitsmass {
  my ($meas, $model) =  @_;
  my $avg=average($meas);
  my $s_meas_2 = sumover(($meas-$avg)**2)/($meas->getdim("0")-1);
  my $s_model2 = sumover(($model-$avg)**2)/($meas->getdim("0")-1);
  my $B = $s_model2/ $s_meas_2;
  return $B;
}

# Bildet die (numerische) Differentiation eine
# Funktion
# Übergeben werden die Funktionswerte in einem Vektor
# Funktion legt (noch) keinen Wert auf Abstände der
# zugehörigen x-Werte, sondern berechnet nur die
# Abstände der Funktionswerte
sub differentiate {
  my ($vector) = @_;
  my $di = $vector->getdim("0")-2;
  my $di2= $vector->getdim("0")-1;
  #print STDOUT  $di . " " . $di2 . "\n";
  my $vector_lower = $vector->copy;
  $vector_lower->slice("1:$di2") .= $vector->slice("0:$di");
  $vector_lower->slice("0") .=0;
  my $d_vector= $vector - $vector_lower;
  return $d_vector;
}

# Funtion für die Anpassung mit "lmfit"
# vgl Dokumentation für PDL::Fit:LM für Parameter
sub alpha_LNVT {
  my ($f,$par,$ym,$dyda) = @_;
  my ($x50,$sigma_ln,$phi) = map {$par->slice("($_)")} (0..2);
  my (@dy) = map {$dyda->slice(",($_)")} (0..2);
  # some constraints on the fit process:
#   if ($x50<0) {$x50=1e-6;}
  # if ($x50>1e-5) {$x50=1e-19;}
#   if ($phi<0) {$phi=0.01;}
#   if ($sigma_ln<0) {$sigma_ln=1;}
#   # $mu und $sigma sind die Werte der Partikelgrößenverteilung
   #print STDOUT ("Before x50 $x50 sigma_ln $sigma_ln phi $phi\n");
  my $dQ =dQ_lognorm_from_parameters_2($xm, $x50, $sigma_ln);
#  print STDOUT ("dQ $dQ\n");
  $ym.=alpha_ges_Faran($f, $alpha_abs_PL1, $alpha_abs_PL2, $alpha_abs_PT1,
           $xm, $dQ, $c_F, $h_L, $h_T, $g,$phi)*$dB
        + $alpha_0 *(1-$phi);
#        + alpha_visc($f,
#                           $xm/2,
#                           $dQ,
#                           $phi, $K_P, $K_F, $rho_P, $rho, $eta_F);
  #if ($x50<0) {$x50.=1e-5;}
  # Ableitungen nach den Parametern
  # dy/dmu
  my $d_x50 = $x50/1000;
  my $x50_a = $x50+$d_x50;
  $dQ =dQ_lognorm_from_parameters_2($xm, $x50_a, $sigma_ln);
  my $result = alpha_ges_Faran($f, $alpha_abs_PL1, $alpha_abs_PL2, $alpha_abs_PT1,
          $xm, $dQ, $c_F, $h_L, $h_T, $g, $phi)*$dB
        + $alpha_0 *(1-$phi);
  #$result += alpha_visc($f,
  #                         $xm/2,
  #                         $dQ,
  #                         $phi, $K_P, $K_F, $rho_P, $rho, $eta_F);
  $x50_a=$x50-$d_x50;
  $dQ =dQ_lognorm_from_parameters_2($xm, $x50_a, $sigma_ln);
  my $result2 = alpha_ges_Faran($f, $alpha_abs_PL1, $alpha_abs_PL2, $alpha_abs_PT1,
          $xm, $dQ, $c_F, $h_L, $h_T, $g, $phi)*$dB
        + $alpha_0 *(1-$phi);
  #$result2 += alpha_visc($f,
  #                         $xm/2,
  #                         $dQ,
  #                         $phi, $K_P, $K_F, $rho_P, $rho, $eta_F);
  $dy[0] .= ($result-$result2)/(2*$d_x50);
  # dy/dsigma
  my $d_sigma_ln = $sigma_ln/1000;
  my $sigma_ln_a=$sigma_ln+$d_sigma_ln;
  $dQ =dQ_lognorm_from_parameters_2($xm, $x50, $sigma_ln_a);
  $result = alpha_ges_Faran($f, $alpha_abs_PL1, $alpha_abs_PL2, $alpha_abs_PT1,
          $xm, $dQ, $c_F, $h_L, $h_T, $g, $phi)*$dB
        + $alpha_0 *(1-$phi);
  #$result += alpha_visc($f,
  #                         $xm/2,
  #                         $dQ,
  #                         $phi, $K_P, $K_F, $rho_P, $rho, $eta_F);
  $sigma_ln_a=$sigma_ln-$d_sigma_ln;
  $dQ =dQ_lognorm_from_parameters_2($xm, $x50, $sigma_ln_a);
  $result2 = alpha_ges_Faran($f, $alpha_abs_PL1, $alpha_abs_PL2, $alpha_abs_PT1,
          $xm, $dQ, $c_F, $h_L, $h_T, $g, $phi)*$dB
        + $alpha_0 *(1-$phi);
  #$result2 += alpha_visc($f,
  #                         $xm/2,
  #                         $dQ,
  #                         $phi, $K_P, $K_F, $rho_P, $rho, $eta_F);
  $dy[1] .= ($result-$result2)/(2*$d_sigma_ln);
  # dy/dc
  my $d_phi = $phi/1000;
  my $phi_a=$phi+$d_phi;
  $dQ =dQ_lognorm_from_parameters_2($xm, $x50, $sigma_ln);
  $result = alpha_ges_Faran($f, $alpha_abs_PL1, $alpha_abs_PL2, $alpha_abs_PT1,
          $xm, $dQ, $c_F, $h_L, $h_T, $g, $phi_a)*$dB
        + $alpha_0 *(1-$phi_a);
  #$result += alpha_visc($f,
  #                         $xm/2,
  #                         $dQ,
  #                         $phi, $K_P, $K_F, $rho_P, $rho, $eta_F);
  $phi_a=$phi-$d_phi;
  $result2 = alpha_ges_Faran($f, $alpha_abs_PL1, $alpha_abs_PL2, $alpha_abs_PT1,
          $xm, $dQ, $c_F, $h_L, $h_T, $g, $phi_a)*$dB
        + $alpha_0 *(1-$phi_a);
  #$result2 += alpha_visc($f,
  #                         $xm/2,
  #                         $dQ,
  #                         $phi, $K_P, $K_F, $rho_P, $rho, $eta_F);
  $dy[2] .= ($result-$result2)/(2*$d_phi);
  print STDOUT ("x50 $x50, sigma_ln $sigma_ln, phi $phi\n");
}

# 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) = @_;
  my $no_Stuetz=100;
  # 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=10**$K * ($fmax**$n);
  #print STDOUT  "x50 " . $x50 . "\n";
  # Bestimme Parameter logarithmischen Reihe
  # (Anpassen, falls Stütstellen-Verteilung zu breit/zu eng)
  # kleinere Argumente = breitere Verteilung
  my $c = log10(20/15);
  # Bestimme x(0) (kleinster Wert der Stützstellenverteilung = Index 0)
  # x50 = x($no_Stuetz/2) (x50 hat den Index no_Stuetz/2)
  my $x_0=$x50/ 10**($no_Stuetz/2*$c);
  # Berechne Verteilung aus den vorangegangenen Werten
  my $x = $x_0 * 10**(sequence($no_Stuetz) * $c) ;
  # Rückgabe. x = Klassenmittelwerte. Rückgabewerte sind in m
  #print STDOUT  "x " . $x . "\n";
  return ($x/1e6 ,$x50/1e6);
}
