package size_distribution;

# Everything related to size distributions

#require 5.000;

use strict;
use warnings;
use vars qw(@ISA);
#use Carp;

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

use lib qw( Module );
use Helos;

#@ISA = qw(Exporter);
#@EXPORT = qw(read_Kext);

sub new {
  my( $this ) = shift;
  my $class = ref($this) || $this;
  my $newInst={};

  $newInst->{XM} = (pdl [0.7 , 1.00, 12e-1, 14e-1, 165e-2, 2.00, 2.40, 2.85, 3.40, 4.00, 4.65, 5.50, 6.75, 8.25, 9.75, 11.50, 13.75, 16.50, 19.50, 23.00, 27.50, 33.00, 39.50, 47.00, 56.00, 67.00, 80.00, 95.00, 113, 135, 161])/1e6;
  $newInst->{XO} = pdl[0];
  $newInst->{DX} = pdl[0];
  $newInst->{DQ} = pdl [0];
  $newInst->{X50} = pdl [0];
  $newInst->{SIGMA_LN} = pdl [0];
  $newInst->{C} = pdl [0];
  $newInst->{Q} = pdl [0];
  bless($newInst, $class);

  $newInst->Initialize;
  return $newInst;
}

sub Initialize {
  my ($this) =@_;
  #printf "Initialize $this\n";
  $this->{XM} = (pdl [0.7 , 1.00, 12e-1, 14e-1, 165e-2, 2.00, 2.40, 2.85, 3.40, 4.00, 4.65, 5.50, 6.75, 8.25, 9.75, 11.50, 13.75, 16.50, 19.50, 23.00, 27.50, 33.00, 39.50, 47.00, 56.00, 67.00, 80.00, 95.00, 113, 135, 161])/1e6;
  $this->{XO} = pdl[0];
  $this->{DX} = pdl[0];
  $this->{DQ} = pdl [0];
  $this->{X50} = pdl [0];
  $this->{SIGMA_LN} = pdl [0];
  $this->{C} = pdl [0];
  $this->{Q} = pdl [0];
}

sub xm {
  my $this = shift;
  if (@_) {
    $this->{XM} = shift;
    $this->calculate_xo;
    $this->calculate_dx_frontend;
  }
  return $this->{XM};
}

sub xm_elem {
  my ($this,$index,$value) = @_;
  if ($value) {
    ${$this->{XM}}->slice("$index") = $value;
    $this->calculate_xo;
    $this->calculate_dx_frontend;
  }
  return $this->{XM}->slice("$index")->sclr;
}

sub calculate_xm {
  my ($this) = @_;
  # Calculate lower class limit
  my $x_lower = $this->{XO}->copy;
  # last element of class vector:
  my $i = $this->{XO}->getdim("0")-1;
  # element before last element of class vector:
  my $i1 = $this->{XO}->getdim("0")-2;
  # fill in data
  $x_lower->slice("0") .= $this->{XO}->slice("0")/2;
  $x_lower->slice("1:$i") .= $this->{XO}->slice("0:$i1");
  $this->{XM}=($this->{XO}+$x_lower)/2;
}

sub xo {
  my $this = shift;
  if (@_) {
    $this->{XO} = shift;
    $this->calculate_xm;
    $this->calculate_dx_frontend;
  }
  return $this->{XO};
}

sub xo_elem {
  my ($this,$index,$value) = @_;
  if ($value) {
    ${$this->{XO}}->slice("$index") = $value;
    $this->calculate_xm;
    $this->calculate_dx_frontend;
  }
  return $this->{XO}->slice("$index")->sclr;
}

sub calculate_xo {
  my ($this) =@_;
  my $di = $this->{XM}->getdim("0")-2;
  my $di2= $this->{XM}->getdim("0")-1;
  $this->{XO} = $this->{XM}->copy;
  $this->{XO}->slice("0") .= ($this->{XM}->slice("1")+$this->{XM}->slice("0"))/2;
  for (my $i=1;$i<=$di2;$i++) {
    my $i_minus_1=$i-1;
    $this->{XO}->slice("$i") .= 2 * $this->{XM}->slice("$i") - $this->{XO}->slice("$i_minus_1");
  }
}

sub dx {
  my $this = shift;
  if (@_) {
    $this->{DX} = shift;
  }
  return $this->{DX};
}

sub dx_elem {
  my ($this,$index,$value) = @_;
  if ($value) {
    ${$this->{DX}}->slice("$index") = $value;}
  return $this->{DX}->slice("$index")->sclr;
}


sub calculate_dx_frontend {
  my ($this) = @_;
  # Calculate dx
  my $di = $this->{XO}->getdim("0")-2;
  my $di2= $this->{XO}->getdim("0")-1;
  my $xo_lower = $this->{XO}->copy;
  $xo_lower->slice("1:$di2").= $this->{XO}->slice("0:$di");
  $xo_lower->slice("0") .=$this->{XO}->slice("0")/2;
  $this->{DX} = $this->{XO} - $xo_lower;
}


sub dQ {
  my $this = shift;
  if (@_) { $this->{DQ} = shift }
  return $this->{DQ};
}

sub dQ_elem {
  my ($this,$index,$value) = @_;
  if ($value) {
    ${$this->{DQ}}->slice("$index") = $value;}
  return $this->{DQ}->slice("$index")->sclr;
}


sub x50 {
  my $this = shift;
  if (@_) { $this->{X50} = shift }
  return $this->{X50};
}

sub sigma_ln {
  my $this = shift;
  if (@_) { $this->{SIGMA_LN} = shift }
  return $this->{SIGMA_LN};
}

sub c {
  my $this = shift;
  if (@_) { $this->{C} = shift }
  return $this->{C};
}

sub dQ_from_params {
   my $this = shift;
   $this->{DQ} = dQ_lognorm_from_parameters_2($this->xm, $this->x50, $this->sigma_ln);
}

sub Q {
  my $this = shift;
  return $this->{DQ}->cumusumover;
}

sub Q_elem {
  my ($this,$index) = @_;
  return $this->{DQ}->cumusumover->slice("$index")->sclr;
}

sub calculate_q {
  my ($this) = @_;
  #$this->dQ_from_params;
  $this->{Q} = $this->{DQ}/$this->{DX};
}

sub q {
  my $this = shift;
  if (@_) {
    $this->{Q} = shift;
  }
  return $this->{Q};
}

sub q_elem {
  my ($this,$index,$value) = @_;
  if ($value) {
    ${$this->{Q}}->slice("$index") = $value;}
  return $this->{Q}->slice("$index")->sclr;
}

# 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 ($this) = @_;
  # Bestimme Parameter logarithmischen Reihe
  # (Anpassen, falls Stütstellen-Verteilung zu breit/zu eng)
  my $c = log10(7/6);
  # Bestimme x(0) (kleinster Wert der Stützstellenverteilung = Index 0)
  # x50 = x(15) (x50 hat den Index 15)
  my $x_0=$this->{X50} / 10**(15*$c);
  # Berechne Verteilung aus den vorangegangenen Werten
  my $x = $x_0 * 10**(sequence(30) * $c) ;
  $this->xm($x);
}



1;
