package Helos;
require 5.000;
require Exporter;
use Carp;

@ISA = qw(Exporter);
@EXPORT = qw(read_Helos calculate_dx xm_from_xo dQ_from_Q dQ2 dQ0 dQ_norm_from_parameters dQ_lognorm_from_parameters dQ_lognorm_from_parameters_2 q3_star);

#!/usr/bin/perl

use strict;
use warnings;


sub read_Helos {
  my ($helos_file)=@_;

  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 $x50;           # PGV: Medianwert in m
  my $x84;           # PGV: x_84-Wert in m
  my $sigma_ln;      # PGV: Standardabweichung der Lognormalverteilung
  my $linebuf;
  my $index=0;
  my @fields;

  # HELOS - Datei einlesen
  open(FILE,$helos_file) or croak("Kann HELOS-Datei $helos_file nicht öffnen \n");
  # Header einlesen, bis der String "dQ3" gefunden wird
  do {
    $linebuf=<FILE>;
    chomp($linebuf);
  } until (index($linebuf,"dQ3")!=-1);
  $linebuf=<FILE>;
  # Daten einlesen und in die entsprechendenen Felder schreiben
  $xo_f[$index]=0;
  while ( (chomp($linebuf=<FILE>) ne "" ) and (length($linebuf)>1)  ){
    @fields=split(" ",$linebuf);
    $xo_f[$index]=$fields[0] / 1e6; # Umrechnung von µm in m
    $Q3_f[$index]=$fields[1] / 100; # Umrechnung von % in Werte von 0-1
    $xm_f[$index]=$fields[3] / 1e6; # Umrechnung von µm in m
    $q3_f[$index]=$fields[4];
    $dQ3_f[$index]=$fields[6] /100; # Umrechnung von % in Werte von 0-1
    $index++;
  }
  #$linebuf=<FILE>;
  # Extract x50
  # value is in µm, convert it to m
  $linebuf=<FILE>;
  chop($linebuf);
  @fields=split(" ",$linebuf);
  my $i=0;
  while ( ($i<=$#fields) && ($fields[$i] ne "x50") ) {$i++;}
  unless ($i>=$#fields) {$x50=$fields[$i+2] / 1e6; }
  # Extract x16, x84
  # values are in µm, convert them to m
  $linebuf=<FILE>;
  chop($linebuf);
  @fields=split(" ",$linebuf);
  $i=0;
  while ( ($i<=$#fields) && ($fields[$i] ne "x84") ) {$i++;}
  unless ($i>=$#fields) {$x84=$fields[$i+2] / 1e6; }
  $sigma_ln = log($x84/$x50);
  close(FILE);
  return (\@xo_f,\@Q3_f,\@xm_f,\@q3_f,\@dQ3_f,$x50,$sigma_ln);
}

# Calculate dx (class size or width) from a given xo
# (upper boundary of classes)

sub calculate_dx {
  use PDL;
  use PDL::Math;
  use PDL::Primitive;
  use PDL::Complex;
  my ($xo) = @_;
  # Calculate dx
  my $di = $xo->getdim("0")-2;
  my $di2= $xo->getdim("0")-1;
  my $xo_lower = $xo->copy;
  $xo_lower->slice("1:$di2").= $xo->slice("0:$di");
  $xo_lower->slice("0") .=$xo->slice("0")/2;
  my $dx = $xo - $xo_lower;
  return $dx;
}

# Calculate xm from a give xo
# (middle of classes from upper class boundaries

sub xm_from_xo {
  use PDL;
  use PDL::Math;
  use PDL::Primitive;
  use PDL::Complex;
  my ($xo)= @_;
  my $di=$xo->getdim("0")-1;
  my $di2=$di-1;
  my $xm=zeroes($di+1);
  $xm->slice("1:$di") .=($xo->slice("0:$di2")->copy+ $xo->slice("1:$di")->copy)/2;
  $xm->slice("0").=2*$xo->slice("0")-$xm->slice("1");
  return($xm);
}


# Calculates the projected-area particle size distribution from
# a given Volume-based distribution
# according to ISO
#
# Parameters: Size distribution dQ3
#             corresponding particle size classes xm (mean of classes)
#
# Returns:    Projected-area-based distribution dQ2 which corresponds to xm

sub dQ2 {
  use PDL;
  use PDL::Math;
  use PDL::Primitive;
  use PDL::Complex;
  use PDL::Graphics2D;
  #use PDL::Graphics::PGPLOT;

  my ($dQ3,$xm)=@_;
  #my $dQ3 = pdl @{$dQ3_f};
  #my $xm = pdl @{$xm_f};
  my $dQ2;         # PGV: dQ2
  $dQ2= ($xm*$dQ3) / sum($xm*$dQ3);
  return $dQ2;
}

# Calculates the number-based particle size distribution from
# a given Volume-based distribution
# according to ISO
#
# Parameters: Size distribution dQ3
#             corresponding particle size classes xm (mean of classes)
#
# Returns:    number-based distribution dQ0 which corresponds to xm


sub dQ0 {
  use PDL;
  use PDL::Math;
  use PDL::Primitive;
  use PDL::Complex;
  use PDL::Graphics2D;
  #use PDL::Graphics::PGPLOT;

  my ($dQ3,$xm)=@_;
  #my $dQ3 = pdl @{$dQ3_f};
  #my $xm = pdl @{$xm_f};
  my $dQ0;         # PGV: dQ2
  $dQ0= ($xm**(-3) * $dQ3) / sum($xm**(-3) * $dQ3);
  return $dQ0;
}

sub dQ_from_Q {
  use PDL;
  use PDL::Math;
  use PDL::Primitive;

  my ($Q) = @_;
  my $Q_shift=$Q->copy->rotate(1);
  $Q_shift->slice("0").=0;
  my $dQ = $Q - $Q_shift;
  return $dQ;
}


# Calculates a normal distribution based upon given parameters
# (mean and standard deviation and number of fields in array)
# returns xm and dQ

sub dQ_norm_from_parameters {
  use PDL;
  use PDL::Math;
  use PDL::Primitive;
  use PDL::Complex;
  use PDL::Graphics2D;
  #use PDL::Graphics::PGPLOT;
  my ($mu, $sigma, $n) = @_;
  # definition of pi
  my $pi = 2 *asin(1);
  # k contains the threshold for limiting the array
  # (borders of array = b1,b2, where q(b1) < k and q(b2)<k)
  my $k = 0.001;
  # Calculate borders
  my $b1 = $mu - $sigma * sqrt(-2*log($k*$sigma*sqrt(2*$pi)));
  my $b2 = $mu + $sigma * sqrt(-2*log($k*$sigma*sqrt(2*$pi)));
  # Calculate x_o (upper limit of classes)
  my $dx = ($b2 - $b1) / $n;
  my $xo = sequence($n+1)* $dx + $b1;
  # Calculate q
  my $q = 1 / ($sigma * sqrt(2*$pi)) * exp(- ($xo-$mu)**2 / (2*$sigma**2));
  # Calculate x_m (center of classes)
  my $xm = $xo->slice("1:$n")-$dx/2;
  # Calculate dQ
  my $dQ = $q->slice("1:$n")*$dx;
  my $Q = sumover($dQ);
  return ($xm,$dQ,$xo, $q);
}

# Calculates a lognormal distribution based upon given parameters
# (mean and standard deviation and number of fields in array)
# returns xm and dQ

sub dQ_lognorm_from_parameters {
  use PDL;
  use PDL::Math;
  use PDL::Primitive;
  use PDL::Complex;
  use PDL::Graphics2D;
  #use PDL::Graphics::PGPLOT;
  my ($mu, $sigma, $n) = @_;
  # definition of pi
  my $pi = 2 *asin(1);
  # k contains the threshold for limiting the array
  # (borders of array = b1,b2, where q(b1) < k and q(b2)<k)
  my $k = 0.01;
  # Calculate borders
  my $b1_ln = - $sigma**2 + log($mu) - sqrt( ($sigma**2-log($mu))**2 + (log($mu))**2 - 2 * $sigma**2 * log($k*$sigma*sqrt(2*$pi)));
  my $b2_ln = - $sigma**2 + log($mu) + sqrt( ($sigma**2-log($mu))**2 + (log($mu))**2 - 2 * $sigma**2 * log($k*$sigma*sqrt(2*$pi)));
  my $b1 = exp($b1_ln);
  my $b2 = exp($b2_ln);
  #my $b2 = $mu + $sigma * sqrt(-2*log($k*$sigma*sqrt(2*$pi)));
  # Calculate x_o (upper limit of classes)
  #my $dx = ($b2 - $b1) / $n;
  #my $xo = sequence($n+1)* $dx + $b1;
  my $dx_ln = ($b2_ln - $b1_ln) / $n;
  my $xo_ln = sequence($n+1)* $dx_ln + $b1_ln;
  my $xo = exp($xo_ln);
  # Calculate q
  my $q = 1 / ($sigma * sqrt(2*$pi)) *(1/$xo)* exp(- (log($xo/$mu))**2 / (2*$sigma**2));
  # Calculate x_m (center of classes)
  #my $xm = $xo->slice("1:$n")-$dx/2;
  my $xo2 = $xo->slice("1:$n");
  my $n_1 = $n-1;
  my $dx = $xo2-$xo->slice("0:$n_1");
  my $xm = $xo->slice("1:$n") - $dx/2;
  # Calculate dQ
  my $dQ = $q->slice("1:$n")*$dx;
  my $Q = sumover($dQ);
  #printf("Q $Q\n");
  return ($xm,$dQ,$xo, $q);
}

# Calculates a lognormal distribution based upon given parameters
# (mean and standard deviation and array of particle sizes)
#

sub dQ_lognorm_from_parameters_2 {
  use PDL;
  use PDL::Math;
  use PDL::Primitive;
  use PDL::Complex;
  use PDL::Graphics2D;
  #use PDL::Graphics::PGPLOT;
  my ($xm, $mu, $sigma) = @_;
  #printf "mu ". $mu . "\n";
  # definition of pi
  my $pi = 2 *asin(1);
  # Calculate q
  my $q = 1 / ($sigma * sqrt(2*$pi)) * (1/$xm) * exp(- (log($xm/$mu))**2 / (2*$sigma**2));
  # Calculate dQ
  my $di = $xm->getdim("0")-2;
  my $di2= $xm->getdim("0")-1;
  my $xo = $xm->copy;
  $xo->slice("0") .= ($xm->slice("1")+$xm->slice("0"))/2;
  for (my $i=1;$i<=$di2;$i++) {
    my $i_min=$i-1;
    $xo->slice("$i") .= 2 * $xm->slice("$i") - $xo->slice("$i_min");
  }
  my $xo_lower = $xo->copy;
  $xo_lower->slice("1:$di2") .= $xo->slice("0:$di");
  $xo_lower->slice("0") .=$xo->slice("0")/2;
  my $dx = $xo - $xo_lower;
  my $dQ = $q*$dx;
  #my $Q = sumover($dQ);
  return ($dQ);
}

# Berechnet die Dichteverteilung des dekadischen Logarithmus (von x)
# Übergeben werden:
# xm: Vektor der Partiklegrößen (Klassenmittelwerte)
# q3: zugehöriger Vektor der Dichteverteilung
# Einheiten: sollten sich aufheben (z.B. x in m, q3 in 1/m)
# Rückgabe Dichteverteilung des dekadischen Logarithmus (von x)
# zur Darstellung im logarithmisch geteilten Diagramm

sub q3_star {
  use PDL;
  use PDL::Math;
  use PDL::Primitive;
  use PDL::Complex;
  my ($xm,$q3)=@_;
  my $q_c_star=$xm/log10(exp(1))*$q3;
  # log10(exp(1)) = lg e
  return $q_c_star;
}

1;
