package Faran_lsg_pdl;
require 5.000;
require Exporter;
use Carp;

@ISA = qw(Exporter);
@EXPORT = qw(k alpha Q_Ext_Faran Q_Ext_Faran_noweight Q_Sca_Faran_noweight alpha_ges_Faran alpha_ges_Faran_noweight alpha_ges_Faran_real j_n j_prime_n y_n L_Faran L_Faran_real F_inf_Faran sigma_s alpha_Str_Faran alpha_Str_Faran_noweight A_Str_Faran c_Faran);

#!/usr/bin/perl

# Faran_lsg_pdl.pm
# (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.


# Faran_lsg_pdl.pm berechnet die Extinktionswirkung eines Partikels
# definierter Größe für gegebene Werte des Brechungsindex (Real-
# und Imaginäranteil). Die Berechnungen erfolgen nach
# Faran 1951, in der von Hay und Mercer (1985) gewählten Notation
#
# Das Modul besitzt eine zentrale Funktion, Q_Ext_Faran
# Dies ist die Funktion, die den Koeffizienten berechnet.
# Der Aufruf der Funktion wird anhand des folgenden Beispiels deutlich:
#
#  use lib qw( modelle); # Pfad zu diesem Modul
#  use Mie_lsg_pdl;
#  use PDL;
#  use PDL::Math;
#  use PDL::Complex;
#  my $frequenzen = pdl [10, 100, 200]; # in MHz
#  my $partikelgroesse = pdl [50/1e6, 148.545/1e6,1/1e6]; # in µm
#  my $res = Q_Ext_Faran($frequenzen,$partikelgroesse,(pdl[1,0,0]),1400,1.33-0.1*i);
#  printf("$res\n");

use strict;
use warnings;

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

my $dB=log(10)/20;

###################################### End of main part ###########
# subroutines following


# Definition des Mie-Parameters

sub alpha {
  my ($x, $f, $c) = @_;
  my $pi=2*asin(1)->sclr;
  #my $alpha = $pi*$x/$c * $f;
  my $alpha = $pi/$c*outer($x,$f);
  return $alpha;
}

# Bessel function
sub j_n {
  my ($ka_param,$n_max, $n) = @_;
  my $ka;
  # test if we have a complex wave number $ka
  if ($ka_param->ndims==2) {
    # No, piddle has 2 dims -> no complex part)
    # convert it to complex number
    #printf "Params " . $ka_param->ndims . "\n";
    #printf "ka " . $ka_param . "\n";
    $ka=r2C($ka_param);
  } else { $ka=$ka_param->copy; }
  # Initialize arrays (dimensions + fill with 0)
  my $j = r2C(zeroes($ka->getdim(1),$ka->getdim(2),$n_max));
  # Calculate first 2 terms
  $j->slice(':,:,:,0').=sin($ka)/$ka;
  $j->slice(':,:,:,1').=sin($ka)/$ka**2 - cos($ka)/$ka;
  # Calculate all other terms based on the preceeding
  # Bessel function series for ka
  for(my $i=2;$i<$n_max;$i++) {
    my $mask = r2C($n>=$i);# + i*($n>=$i);
    #printf "j " . $j->slice(":,:,:,".$i) . "\n";
    #printf "mask " . ($mask*((2*($i-1)+1)/$ka)) . "\n";
    $j->slice(":,:,:,$i").=$mask*((2*($i-1)+1)/$ka*$j->slice(":,:,:,".($i-1))-$j->slice(":,:,:,".($i-2)));
    #printf("i $i --j $j\n");
  }
  return $j;
}

sub j_n_real {
  my ($ka,$n_max, $n) = @_;
  # Initialize arrays (dimensions + fill with 0)
  my $j = zeroes($ka->getdim(0),$ka->getdim(1),$n_max);
  # Calculate first 2 terms
  $j->slice(':,:,0').=sin($ka)/$ka;
  $j->slice(':,:,1').=sin($ka)/$ka**2 - cos($ka)/$ka;
  # Calculate all other terms based on the preceeding
  # Bessel function series for ka
  for(my $i=2;$i<$n_max;$i++) {
    my $mask = $n>=$i;
    # Hack because I don't know how to code expressions (e.g. $i-1)
    # if contained in a slice string
    my $i_minus1= $i-1;
    my $i_minus2= $i-2;
    $j->slice(":,:,$i").=$mask*((2*$i_minus1+1)/$ka*$j->slice(":,:,$i_minus1")-$j->slice(":,:,$i_minus2"));
    #printf("i $i --j $j\n");
  }
  return $j;
}


# Neumann function
sub y_n {
  my ($ka,$n_max, $n) = @_;
  # test if we have a complex wave number $ka
  if ($ka->ndims==2) {
    # No, piddle has 2 dims -> no complex part)
    # convert it to complex number
    $ka=r2C($ka);
  }
  # Initialize arrays (dimensions + fill with 0)
  my $y = r2C(zeroes(double,$ka->getdim(1),$ka->getdim(2),$n_max));
  # Calculate first 2 terms
  $y->slice(':,:,:,0').= - cos($ka)/$ka;
  #$y->slice(':,:,:,0').= 1+1*i;
  $y->slice(':,:,:,1').= - sin($ka)/$ka - cos($ka)/$ka**2;
  #$y->slice(':,:,:,1').= 2+2*i;
  #$ka = 1+1*i;
  # Calculate all other terms based on the preceeding
  for(my $i=2;$i<$n_max;$i++) {
    my $mask = r2C($n>=$i);
    # Hack because I don't know how to code expressions (e.g. $i-1)
    # if contained in a slice string
    my $i_minus1= $i-1;
    my $i_minus2= $i-2;
    $y->slice(":,:,:,$i").=$mask*((2*$i_minus1+1)/$ka*$y->slice(":,:,:,$i_minus1")-$y->slice(":,:,:,$i_minus2"));
    #$y->slice(":,:,:,$i").=((2*$i_minus1+1)/$ka*$y->slice(":,:,:,$i_minus1")-$y->slice(":,:,:,$i_minus2"));
  }
  return $y;
}

sub y_n_real {
  my ($ka,$n_max, $n) = @_;
  # Initialize arrays (dimensions + fill with 0)
  my $y = zeroes(double,$ka->getdim(0),$ka->getdim(1),$n_max);
  # Calculate first 2 terms
  $y->slice(':,:,0').= - cos($ka)/$ka;
  $y->slice(':,:,1').= - sin($ka)/$ka - cos($ka)/$ka**2;
  # Calculate all other terms based on the preceeding
  for(my $i=2;$i<$n_max;$i++) {
    my $mask = $n>=$i;
    # Hack because I don't know how to code expressions (e.g. $i-1)
    # if contained in a slice string
    my $i_minus1= $i-1;
    my $i_minus2= $i-2;
    $y->slice(":,:,$i").=$mask*((2*$i_minus1+1)/$ka*$y->slice(":,:,$i_minus1")-$y->slice(":,:,$i_minus2"));
  }
  return $y;
}


# derivative of the Bessel function
sub j_prime_n {
  my ($ka, $n_max, $n, $j) = @_;
  # test if we have a complex wave number $ka
  if ($ka->ndims==2) {
    # No, piddle has 2 dims -> no complex part)
    # convert it to complex number
    $ka=r2C($ka);
  }
  # Initialize arrays (define dimensions + fill with 0 )
  my $j_prime = r2C(zeroes($ka->getdim(1),$ka->getdim(2),$n_max));
  # Calculate first term
  $j_prime->slice(":,:,:,0").=-$j->slice(':,:,:,1');
  # Calculate remaining terms
  for(my $i=1;$i<$n_max-1;$i++) {
    my $mask = r2C($n>=$i);
    # Hack because I don't know how to code expressions (e.g. $i-1)
    # if contained in a slice string
    my $i_minus= $i-1;
    my $i_plus = $i+1;
    $j_prime->slice(":,:,:,$i").=$mask*
      ( ($i*$j->slice(":,:,:,$i_minus") - ($i+1)*$j->slice(":,:,:,$i_plus")) /
        (2*$i+1) );
  }
  return $j_prime;
}

sub j_prime_n_real {
  my ($ka, $n_max, $n, $j) = @_;
  # Initialize arrays (define dimensions + fill with 0 )
  my $j_prime = zeroes($ka->getdim(0),$ka->getdim(1),$n_max);
  # Calculate first term
  $j_prime->slice(":,:,0").=-$j->slice(':,:,1');
  # Calculate remaining terms
  for(my $i=1;$i<$n_max-1;$i++) {
    my $mask = $n>=$i;
    # Hack because I don't know how to code expressions (e.g. $i-1)
    # if contained in a slice string
    my $i_minus= $i-1;
    my $i_plus = $i+1;
    $j_prime->slice(":,:,$i").=$mask*
      ( ($i*$j->slice(":,:,$i_minus") - ($i+1)*$j->slice(":,:,$i_plus")) /
        (2*$i+1) );
  }
  return $j_prime;
}


# derivative of the Neumann function
sub y_prime_n {
  my ($ka, $n_max, $n, $y) = @_;
  # test if we have a complex wave number $ka
  if ($ka->ndims==2) {
    # No, piddle has 2 dims -> no complex part)
    # convert it to complex number
    $ka=r2C($ka);
  }
  # Initialize arrays (define dimensions + fill with 0 )
  my $y_prime = r2C(zeroes($ka->getdim(1),$ka->getdim(2),$n_max));
  # Calculate first term
  $y_prime->slice(":,:,:,0").=-$y->slice(':,:,:,1');
  # Calculate remaining terms
  for(my $i=1;$i<$n_max-1;$i++) {
    my $mask = r2C($n>=$i);
    # Hack because I don't know how to code expressions (e.g. $i-1)
    # if contained in a slice string
    my $i_minus= $i-1;
    my $i_plus = $i+1;
    $y_prime->slice(":,:,:,$i").=$mask*
      ( ($i*$y->slice(":,:,:,$i_minus") - ($i+1)*$y->slice(":,:,:,$i_plus")) /
        (2*$i+1) );
  }
  return $y_prime;
}

sub y_prime_n_real {
  my ($ka, $n_max, $n, $y) = @_;
  # Initialize arrays (define dimensions + fill with 0 )
  my $y_prime = zeroes($ka->getdim(0),$ka->getdim(1),$n_max);
  # Calculate first term
  $y_prime->slice(":,:,0").=-$y->slice(':,:,1');
  # Calculate remaining terms
  for(my $i=1;$i<$n_max-1;$i++) {
    my $mask = $n>=$i;
    # Hack because I don't know how to code expressions (e.g. $i-1)
    # if contained in a slice string
    my $i_minus= $i-1;
    my $i_plus = $i+1;
    $y_prime->slice(":,:,$i").=$mask*
      ( ($i*$y->slice(":,:,$i_minus") - ($i+1)*$y->slice(":,:,$i_plus")) /
        (2*$i+1) );
  }
  return $y_prime;
}


# Legendre polynom (=solution to Legendre function)
sub P_n {
  my ($theta,$n_max) = @_;
  # Initialize array;
  my $P_n = zeroes($n_max);
  my $mu = cos($theta);
  # calculate first 2 terms
  $P_n->slice("0") .=1;
  $P_n->slice("1") .=cos($theta);
  # remaining terms
  for (my $i=1;$i<$n_max-1;$i++) {
    my $i_plus_1  = $i+1;
    my $i_minus_1 = $i-1;
    $P_n->slice("$i_plus_1") .= ( (2*$i+1)* $mu * $P_n->slice("$i") - $i* $P_n->slice("$i_minus_1") ) /
                                ( $i+1 );
  }
  return $P_n;
}

# prints the first four dimensions of a piddle

sub pdims {
  my ($piddle, $name) = @_;
  my $a=$piddle->getdim("0");
  printf("$name $a ");
  $a=$piddle->getdim("1");
  printf("$a ");
  $a=$piddle->getdim("2");
  printf("$a ");
  $a=$piddle->getdim("3");
  printf("$a\n");
}

# Calculation of the "complex phase shift" (Hay & Mercer)
# Notation adapted from Feuillade
# Feuillade    Hay & Mercer
# ka      x
# ka_PL      x'
# ka_PT      s'

sub tan_eta_n {
  my ($ka, $ka_PL, $ka_PT, $n, $n_PL, $n_PT, $n_max, $g) =@_;
  # Calculate n series
  my $n_ser = sequence($n_max);
  $n_ser = $n_ser->dummy(0,$ka->getdim(2));
  $n_ser = $n_ser->dummy(0,$ka->getdim(1));
  #printf("n_ser $n_ser\n");
  # Calculate Bessel and Neumann function series terms
  my $j = j_n($ka,$n_max,$n);
  my $y = y_n($ka,$n_max,$n);
  my $j_prime = j_prime_n($ka,$n_max,$n,$j);
  my $y_prime = y_prime_n($ka,$n_max,$n,$y);
  # Calculation of tan_delta_n
  my $tan_delta_n = -$j / $y;
  #pdims($tan_delta_n,"tan_delta_n");
  # Calculation of tan_alpha_n
  my $tan_alpha_n = -$ka * $j_prime / $j;
  #pdims($tan_alpha_n,"tan_alpha_n");
  # Calculation of tan_beta_n
  my $tan_beta_n = -$ka * $y_prime / $y;
  #pdims($tan_beta_n,"tan_beta_n");
  # Calculation of material properties coefficient Phi_n
  # tan_alpha_n_PL
  my $j_PL = j_n($ka_PL,$n_max,$n_PL);
  my $j_prime_PL = j_prime_n($ka_PL,$n_max,$n_PL,$j_PL);
  my $tan_alpha_n_PL = - $ka_PL * $j_prime_PL / $j_PL;
  #pdims($tan_alpha_n_PL,"tan_alpha_n_PL");
  # tan_alpha_n_PT
  my $j_PT = j_n($ka_PT,$n_max,$n_PT);
  my $j_prime_PT = j_prime_n($ka_PT,$n_max,$n_PT,$j_PT);
  my $tan_alpha_n_PT = - $ka_PT * $j_prime_PT / $j_PT;
  #pdims($tan_alpha_n_PT,"tan_alpha_n_PT");
  # Phi_n
  #pdims($nom,"nom");
  #pdims($denom, "denom");
  my $Phi_n = ($ka_PT**2) /2 /$g * (
        ($tan_alpha_n_PL / ($tan_alpha_n_PL + 1) -
          (($n_ser**2) + $n_ser)/($n_ser**2 + $n_ser - 1 - ($ka_PT**2)/2 + $tan_alpha_n_PT)) /
        ( (($n_ser**2) + $n_ser  - ($ka_PT**2)/2 + 2 * $tan_alpha_n_PL)/( $tan_alpha_n_PL + 1)
         -( (($n_ser**2) + $n_ser) * ($tan_alpha_n_PT + 1) ) /
            (($n_ser**2) + $n_ser - 1 - ($ka_PT**2)/2 + $tan_alpha_n_PT) )
        );
  #pdims($Phi_n,"Phi_n");
  # eta_n
  my $taneta_n = $tan_delta_n * ($tan_alpha_n + $Phi_n) / ($tan_beta_n + $Phi_n);
  #pdims($taneta_n,"taneta_n");
  $taneta_n->inplace->badmask(0);
  #printf $taneta_n->slice(":,:,:," . ($n_max-10)) . "\n";
#   printf "j_n".  $j->slice(":,:,:,:") . "\n";
#   printf "y_n".  $y->slice(":,:,:,:") . "\n";
#   printf "j_prime_n".  $j_prime->slice(":,:,:,:") . "\n";
#   printf "y_prime_n".  $y_prime->slice(":,:,:,:") . "\n";
#   printf "j_PL_n" . $j_PL->slice(":,:,:,:") . "\n";
#   printf "j_prime_PL_n".  $j_prime_PL->slice(":,:,:,:") . "\n";
#   printf "j_PT_n" . $j_PT->slice(":,:,:,:") . "\n";
#   printf "j_prime_PT_n".  $j_prime_PT->slice(":,:,:,:") . "\n";
  return $taneta_n;
}

sub tan_eta_n_real {
  my ($ka, $ka_PL, $ka_PT, $n, $n_PL, $n_PT, $n_max, $g) =@_;
  # Calculate n series
  my $n_ser = sequence($n_max);
  $n_ser = $n_ser->dummy(0,$ka->getdim(1));
  $n_ser = $n_ser->dummy(0,$ka->getdim(0));
  #printf("n_ser $n_ser\n");
  # Calculate Bessel and Neumann function series terms
  my $j = j_n_real($ka,$n_max,$n);
  my $y = y_n_real($ka,$n_max,$n);
  my $j_prime = j_prime_n_real($ka,$n_max,$n,$j);
  my $y_prime = y_prime_n_real($ka,$n_max,$n,$y);
  # Calculation of tan_delta_n
  my $tan_delta_n = -$j / $y;
  #pdims($tan_delta_n,"tan_delta_n");
  # Calculation of tan_alpha_n
  my $tan_alpha_n = -$ka * $j_prime / $j;
  #pdims($tan_alpha_n,"tan_alpha_n");
  # Calculation of tan_beta_n
  my $tan_beta_n = -$ka * $y_prime / $y;
  #pdims($tan_beta_n,"tan_beta_n");
  # Calculation of material properties coefficient Phi_n
  # tan_alpha_n_PL
  my $j_PL = j_n_real($ka_PL,$n_max,$n_PL);
  my $j_prime_PL = j_prime_n_real($ka_PL,$n_max,$n_PL,$j_PL);
  my $tan_alpha_n_PL = - $ka_PL * $j_prime_PL / $j_PL;
  #pdims($tan_alpha_n_PL,"tan_alpha_n_PL");
  # tan_alpha_n_PT
  my $j_PT = j_n_real($ka_PT,$n_max,$n_PT);
  my $j_prime_PT = j_prime_n_real($ka_PT,$n_max,$n_PT,$j_PT);
  my $tan_alpha_n_PT = - $ka_PT * $j_prime_PT / $j_PT;
  #pdims($tan_alpha_n_PT,"tan_alpha_n_PT");
  # Phi_n
  #pdims($nom,"nom");
  #pdims($denom, "denom");
  my $Phi_n = ($ka_PT**2) /2 /$g * (
        ($tan_alpha_n_PL / ($tan_alpha_n_PL + 1) -
          (($n_ser**2) + $n_ser)/($n_ser**2 + $n_ser - 1 - ($ka_PT**2)/2 + $tan_alpha_n_PT)) /
        ( (($n_ser**2) + $n_ser  - ($ka_PT**2)/2 + 2 * $tan_alpha_n_PL)/( $tan_alpha_n_PL + 1)
         -( (($n_ser**2) + $n_ser) * ($tan_alpha_n_PT + 1) ) /
            (($n_ser**2) + $n_ser - 1 - ($ka_PT**2)/2 + $tan_alpha_n_PT) )
        );
  #pdims($Phi_n,"Phi_n");
  # eta_n
  my $taneta_n = $tan_delta_n * ($tan_alpha_n + $Phi_n) / ($tan_beta_n + $Phi_n);
  #pdims($taneta_n,"taneta_n");
  $taneta_n->inplace->badmask(0);
  #printf $taneta_n->slice(":,:," . ($n_max-10)) . "\n";
  return $taneta_n;
}


sub A_n {
  my ($ka, $ka_PL, $ka_PT, $n, $n_PL, $n_PT, $n_max, $g) =@_;
  my $tan_eta_n = tan_eta_n($ka, $ka_PL, $ka_PT, $n, $n_PL, $n_PT, $n_max, $g);
  my $An = - $tan_eta_n * i / (1+ $tan_eta_n * i);
  #printf "A_n[1] ".  $An->slice(":,:,:,:") . "\n";
  return $An;
}

sub A_n_real {
  my ($ka, $ka_PL, $ka_PT, $n, $n_PL, $n_PT, $n_max, $g) =@_;
  my $tan_eta_n = tan_eta_n_real($ka, $ka_PL, $ka_PT, $n, $n_PL, $n_PT, $n_max, $g);
  my $An = - $tan_eta_n * i / (1+ $tan_eta_n * i);
  #printf $An->slice(":,:,:,1") . "\n";
  return $An;
}



# series termination according to Wiscombe

sub n_B {
  my ($alpha)= @_;
  #my $res = ceil($alpha+ 4.05 * ($alpha ** (1/3)) + 2); # für 8  < alpha < 4200
  #my $res = ceil(1+0.7*($alpha)**0.95)+10; # von Frank, + Summand 10
  # ceil($akpha + 4*($alpha ** (1/3))+1); # für 0.02 <= $alpha <= 8
  # ceil($akpha + 4*($alpha ** (1/3))+2); # für 4200 <= $alpha <= 20000
  #my $res = ceil($alpha)+500;
  my $res= abs(ceil(re($alpha)/re($alpha)*30));
  #my $res = ceil($alpha)/ceil($alpha)*80;
  return $res;
}

# Calculation of Extinction efficiencies

sub Q_Ext_Faran {
 #my ($k, $m_r)= @_;
  my ($f, $xm, $dQ2, $c, $h_L, $h_T, $g)=@_;
  # Calculate the wave numbers
  my $pi = 2*asin(1);
  my $ka=outer($xm,$f)*1e6/$c*$pi; # matmult ?
  my $ka_PL = $ka/$h_L;
  my $ka_PT = $ka/$h_T;
  #Calculate required length of Bessel function series
  my $n_b1 = n_B($ka);
  my $n_max = max($n_b1);
  my $n_b2 = n_B($ka_PL);
  my $n_max2 = max($n_b2);
  if ($n_max2 > $n_max) { $n_max = $n_max2;}
  my $n_b3 = n_B($ka_PT);
  my $n_max3 = max($n_b3);
  if ($n_max3 > $n_max) { $n_max = $n_max3;}
  # Calculate A_n (Kräuter (Diss), Hay & Mercer (1975) )
  my $An = A_n($ka ,$ka_PL ,$ka_PT ,$n_b1 ,$n_b2 ,$n_b3 ,$n_max ,$g);
  # Create Sum
  # 2*n+1
  my $n=sequence($n_max);
  $n = 2*$n+1;
  $n = $n->dummy(0,$ka->getdim(1));
  $n = $n->dummy(0,$ka->getdim(0));
  # (2n+1)*re(R_n)
  my $Rn2 = re($An)*$n;
  my $maxdim=$n_max-1;
  $Rn2 = $Rn2->slice(":,:,0:$maxdim");
  # Sum the above
  my $sum=dsumover($Rn2->xchg(2,0));
  # make sure the dimensions are in the right order
  $sum = $sum->xchg(0,1);
  # divide by ka**2, multiply by 2
  my $Kext_array = 4*$sum/($ka**2);
  # multiply with particle size distribution
  my $weighted = -inner($Kext_array,$dQ2);
  #printf("weighted $weighted\n");
  return $weighted;
}

# same as above, but without weighting
sub Q_Ext_Faran_noweight {
 #my ($k, $m_r)= @_;
  my ($f, $xm, $c, $h_L, $h_T, $g)=@_;
  # Calculate the wave numbers
  my $pi = 2*asin(1);
  my $ka=outer($xm,$f)*1e6/$c*$pi; # matmult ?
  my $ka_PL = $ka/$h_L;
  my $ka_PT = $ka/$h_T;
  #Calculate required length of Bessel function series
  my $n_b1 = n_B($ka);
  my $n_max = max($n_b1);
  my $n_b2 = n_B($ka_PL);
  my $n_max2 = max($n_b2);
  if ($n_max2 > $n_max) { $n_max = $n_max2;}
  my $n_b3 = n_B($ka_PT);
  my $n_max3 = max($n_b3);
  if ($n_max3 > $n_max) { $n_max = $n_max3;}
  # Calculate A_n (Kräuter (Diss), Hay & Mercer (1975) )
  my $An = A_n($ka ,$ka_PL ,$ka_PT ,$n_b1 ,$n_b2 ,$n_b3 ,$n_max ,$g);
  # Create Sum
  # 2*n+1
  my $n=sequence($n_max);
  $n = 2*$n+1;
  $n = $n->dummy(0,$ka->getdim(1));
  $n = $n->dummy(0,$ka->getdim(0));
  # (2n+1)*re(R_n)
  my $Rn2 = re($An)*$n;
  my $maxdim=$n_max-1;
  $Rn2 = $Rn2->slice(":,:,0:$maxdim");
  # Sum the above
  my $sum=dsumover($Rn2->xchg(2,0));
  # make sure the dimensions are in the right order
  $sum = $sum->xchg(0,1);
  # divide by ka**2, multiply by 2
  my $Kext_array = 4*$sum/($ka**2);
  return $Kext_array;
}


# Scattering coefficient without weighting
sub Q_Sca_Faran_noweight {
 #my ($k, $m_r)= @_;
  my ($f, $xm, $c, $h_L, $h_T, $g)=@_;
  # Calculate the wave numbers
  my $pi = 2*asin(1);
  my $ka=r2C(outer($xm,$f)*1e6/$c*$pi); # matmult ?
  my $ka_PL = $ka/$h_L;
  my $ka_PT = $ka/$h_T;
  #Calculate required length of Bessel function series
  my $n_b1 = n_B($ka);
  my $n_max = max($n_b1);
  my $n_b2 = n_B($ka_PL);
  my $n_max2 = max($n_b2);
  if ($n_max2 > $n_max) { $n_max = $n_max2;}
  my $n_b3 = n_B($ka_PT);
  my $n_max3 = max($n_b3);
  if ($n_max3 > $n_max) { $n_max = $n_max3;}
  # Calculate A_n (Kräuter (Diss), Hay & Mercer (1975) )
  my $An = A_n($ka ,$ka_PL ,$ka_PT ,$n_b1 ,$n_b2 ,$n_b3 ,$n_max ,$g);
  # Create Sum
  # 2*n+1
  my $n=sequence($n_max);
  $n = 2*$n+1;
  $n = $n->dummy(0,$ka->getdim(2));
  $n = $n->dummy(0,$ka->getdim(1));
  # (2n+1)*re(R_n)
  my $Rn2 = (re($An)**2 + im($An)**2)*$n;
  my $maxdim=$n_max-1;
  $Rn2 = $Rn2->slice(":,:,0:$maxdim");
  # Sum the above
  my $sum=dsumover($Rn2->xchg(2,0));
  # make sure the dimensions are in the right order
  $sum = $sum->xchg(0,1);
  # divide by ka**2, multiply by 2
  my $Ksca_array = 4*$sum/($ka**2);
  return $Ksca_array;
}


# calculates the (complex) dispersion wave number based
# on the Waterman & Truell approach (forward and backward scattering
# amplitude

sub k {
  my ($f, $alpha_abs, $xm, $dQ, $c, $h_L, $h_T, $g, $phi)=@_;
  # Calculate the wave numbers
  my $pi = 2*asin(1);
  my $ka=outer($xm,$f)*1e6/$c*$pi #;#  matmult ?
   - i*outer($xm,$alpha_abs)/2;
  my $k = 2*$pi*$f*1e6/$c
    - i*$alpha_abs;
  my $F_0  =L_Faran($f, $alpha_abs, $xm, $c, $h_L, $h_T, $g,0);
  my $F_pi =L_Faran($f, $alpha_abs, $xm, $c, $h_L, $h_T, $g,$pi);
#   printf "F_0  " . $F_0->getdim(0) . "\n";
#   printf $F_0->getdim(1) . "\n";
#   printf "\n";
#   printf "ka  " . $ka->getdim(0) . "\n";
#   printf $ka->getdim(1) . "\n";
#   printf "xm  " . $xm->getdim(0) . "\n";
#   printf $xm->getdim(1) . "\n";
#   printf "dQ " . $dQ->getdim(0) . "\n";
#  my $dummy_0=$F_0/re($ka)/$xm;
#  my $dummy_pi=$F_0/re($ka)/$xm;
  my $dummy_0=$F_0/$ka/r2C($xm);
  my $dummy_pi=$F_0/$ka/r2C($xm);
#   printf "F_0/ka/r2C(xm) " . $dummy_0->getdim(0) . "\n";
#   printf $dummy_0->getdim(1) . "\n";
#   printf $dummy_0->getdim(2) . "\n";
#   printf $dummy_0->getdim(3) . "\n";
  my $sum_0 =inner($dummy_0->transpose,$dQ);
  my $sum_pi=inner($dummy_pi->transpose,$dQ);
  my $kappa=$k*sqrt( (1+3*$phi*$sum_0)**2 - (3*$phi*$sum_pi)**2  );
  return $kappa;
}

# Function based on Feuillade
# calculates the Acoustic scattering length
# equal to "scattering amplitude" (Mobley et al, 1999)

sub L_Faran {
  my ($f, $alpha_abs, $xm, $c, $h_L, $h_T, $g, $theta)=@_;
  # Calculate the wave numbers
  my $pi = 2*asin(1);
  my $ka=outer($xm,$f)*1e6/$c*$pi #; matmult ?
   - i*outer($xm,$alpha_abs)/2;
  #my $ka_PL = $ka/$h_L;
  my $c_PL=$h_L*$c;
  my $ka_PL=outer($xm,$f)*1e6/$c_PL*$pi
   - i*outer($xm,$alpha_abs)/2;
  #my $ka_PT = $ka/$h_T;
  my $c_PT=$h_T*$c;
  my $ka_PT=outer($xm,$f)*1e6/$c_PT*$pi
   - i*outer($xm,$alpha_abs)/2;
  #Calculate required length of Bessel function series
  my $n_b1 = n_B($ka);
  my $n_max = max($n_b1);
  my $n_b2 = n_B($ka_PL);
  my $n_max2 = max($n_b2);
  if ($n_max2 > $n_max) { $n_max = $n_max2;}
  my $n_b3 = n_B($ka_PT);
  my $n_max3 = max($n_b3);
  #if ($n_max3 > $n_max) { $n_max = $n_max3;}
  # Calculate A_n (Kräuter (Diss), Hay & Mercer (1975) )
  my $An = A_n($ka ,$ka_PL ,$ka_PT ,$n_b1 ,$n_b2 ,$n_b3 ,$n_max ,$g);
  # Calculate Legendre Polynomial
  my $P_n = P_n($theta,$n_max);
  $P_n = $P_n->dummy(0,$ka->getdim(2));
  $P_n = $P_n->dummy(0,$ka->getdim(1));
  # Create Sum
  # 2*n+1
  my $n=sequence($n_max);
  $n = (2*$n)+1;
  $n = $n->dummy(0,$ka->getdim(2));
  $n = $n->dummy(0,$ka->getdim(1));
  #printf("n $n\n");
  # (2n+1)*re(R_n)*P_n
  my $Rn2 = $An * $n * $P_n;
  my $maxdim=$n_max-1;
  $Rn2 = $Rn2->slice(":,:,:,0:$maxdim");
  # Sum the above
  my $sum=sumover($Rn2->xchg(3,0));
  # make sure the dimensions are in the right order
  $sum = $sum->xchg(0,1);
  $sum = $sum->xchg(0,2);
  # divide by ka, multiply by a, multiply by i
  my $L = $sum/($ka)*($xm/2)*i;
#   my $L_abs = sqrt( (re($L))**2 + (im($L))**2 );
#   return $L_abs;
  return $L;
}

sub L_Faran_real {
  my ($f, $xm, $c, $h_L, $h_T, $g, $theta)=@_;
  # Calculate the wave numbers
  my $pi = 2*asin(1);
  my $ka=outer($xm,$f)*1e6/$c*$pi; #  matmult ?
  #my $ka_PL = $ka/$h_L;
  my $c_PL=$h_L*$c;
  my $ka_PL=outer($xm,$f)*1e6/$c_PL*$pi;
  #my $ka_PT = $ka/$h_T;
  my $c_PT=$h_T*$c;
  my $ka_PT=outer($xm,$f)*1e6/$c_PT*$pi;
  #Calculate required length of Bessel function series
  my $n_b1 = n_B($ka);
  my $n_max = max($n_b1);
  my $n_b2 = n_B($ka_PL);
  my $n_max2 = max($n_b2);
  if ($n_max2 > $n_max) { $n_max = $n_max2;}
  my $n_b3 = n_B($ka_PT);
  my $n_max3 = max($n_b3);
  #if ($n_max3 > $n_max) { $n_max = $n_max3;}
  # Calculate A_n (Kräuter (Diss), Hay & Mercer (1975) )
  my $An = A_n_real($ka ,$ka_PL ,$ka_PT ,$n_b1 ,$n_b2 ,$n_b3 ,$n_max ,$g);
  # Calculate Legendre Polynomial
  my $P_n = P_n($theta,$n_max);
  $P_n = $P_n->dummy(0,$ka->getdim(1));
  $P_n = $P_n->dummy(0,$ka->getdim(0));
  # Create Sum
  # 2*n+1
  my $n=sequence($n_max);
  $n = (2*$n)+1;
  $n = $n->dummy(0,$ka->getdim(1));
  $n = $n->dummy(0,$ka->getdim(0));
  #printf("n $n\n");
  # (2n+1)*re(R_n)*P_n
  my $Rn2 = $An * $n * $P_n;
  my $maxdim=$n_max-1;
  $Rn2 = $Rn2->slice(":,:,:,0:$maxdim");
  # Sum the above
  my $sum=sumover($Rn2->xchg(3,0));
  # make sure the dimensions are in the right order
  $sum = $sum->xchg(0,1);
  $sum = $sum->xchg(0,2);
  # divide by ka**2, multiply by a, multiply by i
  my $L = $sum/($ka)*($xm/2)*i;
  my $L_abs = sqrt( (re($L))**2 + (im($L))**2 );
  return $L_abs;
}


# "Form function" according to Feuillade, 1999

sub F_inf_Faran {
  my ($f, $xm, $dQ2, $c, $h_L, $h_T, $g, $theta)=@_;
  my $L = L_Faran($f,$xm,$dQ2,$c,$h_L,$h_T,$g,$theta);
  printf("L $L\n");
  my $F = $L / ($xm/4);
  # multiply with particle size distribution
  my $weighted = inner($F,$dQ2);
  #printf("weighted $weighted\n");
  return $weighted;
}

# Extinction cross section according to Feuillade, 1999
sub sigma_s {
  my ($f,$xm,$dQ2,$c,$h_L, $h_T,$g,$theta)=@_;
  my $abs_L = L_Faran($f,$xm,$dQ2,$c,$h_L,$h_T,$g,$theta);
  my $sigma = $abs_L**2;
  # multiply with particle size distribution
  my $weighted = inner($sigma,$dQ2);
  #printf("weighted $weighted\n");
  return $weighted;
}

# Returns the attenuation in Np/m, based upon Parameters
# xm : Particle diameters, in m
# c : speed of sound in the fluid in m/s

sub alpha_ges_Faran {
  my ($f, $alpha_abs_F, $alpha_abs_PL, $alpha_abs_PT, $xm, $dQ2, $c, $h_L, $h_T, $g, $phi)=@_;
  # Calculate the wave numbers
  my $pi = 2*asin(1);
  # fluid
  my $ka=outer($xm,$f)*1e6/$c*$pi #; # matmult ?
   - i*outer($xm,$alpha_abs_F)/2;
  #printf "ka " . $ka . "\n";
  #my $ka_PL = $ka/$h_L;
  my $c_PL=$h_L*$c;
  # particle (compressional/longitudinal) wave
  #printf "c_PL   " . $c_PL . "\n";
  my $ka_PL=outer($xm,$f)*1e6/$c_PL*$pi
   - i*outer($xm,$alpha_abs_PL)/2;
  #my $ka_PT = $ka/$h_T;
  my $c_PT=$h_T*$c;
  # particle (shear/transversal) wave
  my $ka_PT=outer($xm,$f)*1e6/$c_PT*$pi
   - i*outer($xm,$alpha_abs_PT)/2;
  #printf "c_PT   " . $c_PT . "\n";
  #printf "a_PT   " . $alpha_abs_PT . "\n";
  #Calculate required length of Bessel function series
  my $n_b1 = n_B($ka);
  my $n_max = max($n_b1);
  my $n_b2 = n_B($ka_PL);
  my $n_max2 = max($n_b2);
  if ($n_max2 > $n_max) { $n_max = $n_max2;}
  my $n_b3 = n_B($ka_PT);
  #printf "$ka_PL\n";
  my $n_max3 = max($n_b3);
  if ($n_max3 > $n_max) { $n_max = $n_max3;}
  # Calculate A_n (Kräuter (Diss), Hay & Mercer (1975) )
  my $An = A_n($ka ,$ka_PL ,$ka_PT ,$n_b1 ,$n_b2 ,$n_b3 ,$n_max ,$g);
  # Create Sum
  # 2*n+1
  my $n=sequence($n_max);
  $n = 2*$n+1;
  $n = $n->dummy(0,$ka->getdim(2));
  $n = $n->dummy(0,$ka->getdim(1));
  # (2n+1)*re(R_n)
  my $Rn2 = re($An)*$n;
  my $maxdim=$n_max-1;
  $Rn2 = $Rn2->slice(":,:,0:$maxdim");
  # Sum the above
  my $sum=dsumover($Rn2->xchg(2,0));
  # make sure the dimensions are in the right order
  $sum = $sum->xchg(0,1);
  # create alpha
  my $alpha = 3*$phi*$sum/(re($ka)**2);
  #printf re($ka) . "\n";
  # multiply with particle size distribution
  my $weighted = -inner($alpha/$xm,$dQ2);
  #printf("weighted $weighted\n");
  # Divide by Frequency
  my $alpha_div_f = $weighted/ $f;
  # printf("div f $alpha_div_f\n");
  # Result is in Np/m/MHz
  return $alpha_div_f;
}

sub alpha_ges_Faran_real {
  my ($f, $xm, $dQ2, $c, $h_L, $h_T, $g, $phi)=@_;
  # Calculate the wave numbers
  my $pi = 2*asin(1);
  my $ka=outer($xm,$f)*1e6/$c*$pi; # matmult ?
  # + i* outer($xm,$alpha_abs)/2;
  #my $ka_PL = $ka/$h_L;
  my $c_PL=$h_L*$c;
  my $ka_PL=outer($xm,$f)*1e6/$c_PL*$pi;
  # + i* outer($xm,$alpha_abs)/2;
  #my $ka_PT = $ka/$h_T;
  my $c_PT=$h_T*$c;
  my $ka_PT=outer($xm,$f)*1e6/$c_PT*$pi;
  # + i* outer($xm,$alpha_abs)/2;
  #Calculate required length of Bessel function series
  my $n_b1 = n_B($ka);
  my $n_max = max($n_b1);
  my $n_b2 = n_B($ka_PL);
  my $n_max2 = max($n_b2);
  if ($n_max2 > $n_max) { $n_max = $n_max2;}
  my $n_b3 = n_B($ka_PT);
  my $n_max3 = max($n_b3);
  if ($n_max3 > $n_max) { $n_max = $n_max3;}
  # Calculate A_n (Kräuter (Diss), Hay & Mercer (1975) )
  my $An = A_n_real($ka ,$ka_PL ,$ka_PT ,$n_b1 ,$n_b2 ,$n_b3 ,$n_max ,$g);
  # Create Sum
  # 2*n+1
  my $n=sequence($n_max);
  $n = 2*$n+1;
  $n = $n->dummy(0,$ka->getdim(1));
  $n = $n->dummy(0,$ka->getdim(0));
  # (2n+1)*re(R_n)
  my $Rn2 = re($An)*$n;
  my $maxdim=$n_max-1;
  $Rn2 = $Rn2->slice(":,:,0:$maxdim");
  # Sum the above
  my $sum=dsumover($Rn2->xchg(2,0));
  # make sure the dimensions are in the right order
  $sum = $sum->xchg(0,1);
  # create alpha
  my $alpha = 3*$phi*$sum/($ka**2);
  #printf $ka . "\n";
  # multiply with particle size distribution
  my $weighted = -inner($alpha/$xm,$dQ2);
  #printf("weighted $weighted\n");
  # Divide by Frequency
  my $alpha_div_f = $weighted/ $f;
  # printf("div f $alpha_div_f\n");
  # Result is in Np/m/MHz
  return $alpha_div_f;
}

sub alpha_ges_Faran_noweight {
  my ($f, $alpha_abs_F, $alpha_abs_PL, $alpha_abs_PT, $xm, $c, $h_L, $h_T, $g, $phi)=@_;
  # Calculate the wave numbers
  my $pi = 2*asin(1);
  # fluid
  my $ka=outer($xm,$f)*1e6/$c*$pi #; # matmult ?
   - i*outer($xm,$alpha_abs_F)/2;
  #printf "ka " . $ka . "\n";
  #my $ka_PL = $ka/$h_L;
  my $c_PL=$h_L*$c;
  # particle (compressional/longitudinal) wave
  #printf "c_PL   " . $c_PL . "\n";
  my $ka_PL=outer($xm,$f)*1e6/$c_PL*$pi
   - i*outer($xm,$alpha_abs_PL)/2;
  #my $ka_PT = $ka/$h_T;
  my $c_PT=$h_T*$c;
  # particle (shear/transversal) wave
  my $ka_PT=outer($xm,$f)*1e6/$c_PT*$pi
   - i*outer($xm,$alpha_abs_PT)/2;
  #printf "c_PT   " . $c_PT . "\n";
  #printf "a_PT   " . $alpha_abs_PT . "\n";
  #Calculate required length of Bessel function series
  my $n_b1 = n_B($ka);
  my $n_max = max($n_b1);
  my $n_b2 = n_B($ka_PL);
  my $n_max2 = max($n_b2);
  if ($n_max2 > $n_max) { $n_max = $n_max2;}
  my $n_b3 = n_B($ka_PT);
  #printf "$ka_PL\n";
  my $n_max3 = max($n_b3);
  if ($n_max3 > $n_max) { $n_max = $n_max3;}
  # Calculate A_n (Kräuter (Diss), Hay & Mercer (1975) )
  my $An = A_n($ka ,$ka_PL ,$ka_PT ,$n_b1 ,$n_b2 ,$n_b3 ,$n_max ,$g);
  # Create Sum
  # 2*n+1
  my $n=sequence($n_max);
  $n = 2*$n+1;
  $n = $n->dummy(0,$ka->getdim(2));
  $n = $n->dummy(0,$ka->getdim(1));
  # (2n+1)*re(R_n)
  my $Rn2 = re($An)*$n;
  my $maxdim=$n_max-1;
  $Rn2 = $Rn2->slice(":,:,0:$maxdim");
  # Sum the above
  my $sum=dsumover($Rn2->xchg(2,0));
  # make sure the dimensions are in the right order
  $sum = $sum->xchg(0,1);
  # create alpha
  my $alpha = -3*$phi*$sum/(re($ka)**2);
  #printf re($ka) . "\n";
  # Divide by Frequency
  my $alpha_div_f =transpose(transpose($alpha/$xm)/$f);
  # printf("div f $alpha_div_f\n");
  # Result is in Np/m/MHz
  return $alpha_div_f;
}


sub c_Faran {
  my ($f, $alpha_abs, $xm, $dQ2, $c, $h_L, $h_T, $g, $phi)=@_;
  # Calculate the wave numbers
  my $pi = 2*asin(1);
  my $omega= 2*$pi*$f*1e6;
  my $ka=outer($xm,$f)*1e6/$c*$pi#; # matmult ?
   - i*outer($xm,$alpha_abs)/2;
  #my $ka_PL = $ka/$h_L;
  my $c_PL=$h_L*$c;
  my $ka_PL=outer($xm,$f)*1e6/$c_PL*$pi
   - i*outer($xm,$alpha_abs)/2;
  #my $ka_PT = $ka/$h_T;
  my $c_PT=$h_T*$c;
  my $ka_PT=outer($xm,$f)*1e6/$c_PT*$pi
   - i*outer($xm,$alpha_abs)/2;
  #Calculate required length of Bessel function series
  my $n_b1 = n_B($ka);
  my $n_max = max($n_b1);
  my $n_b2 = n_B($ka_PL);
  my $n_max2 = max($n_b2);
  if ($n_max2 > $n_max) { $n_max = $n_max2;}
  my $n_b3 = n_B($ka_PT);
  #printf "$ka_PL\n";
  my $n_max3 = max($n_b3);
  if ($n_max3 > $n_max) { $n_max = $n_max3;}
  # Calculate A_n (Kräuter (Diss), Hay & Mercer (1975) )
  my $An = A_n($ka ,$ka_PL ,$ka_PT ,$n_b1 ,$n_b2 ,$n_b3 ,$n_max ,$g);
  # Create Sum
  # 2*n+1
  my $n=sequence($n_max);
  $n = 2*$n+1;
  $n = $n->dummy(0,$ka->getdim(2));
  $n = $n->dummy(0,$ka->getdim(1));
  # (2n+1)*re(R_n)
  my $Rn2 = im($An)*$n;
  my $maxdim=$n_max-1;
  $Rn2 = $Rn2->slice(":,:,0:$maxdim");
  # Sum the above
  my $sum=dsumover($Rn2->xchg(2,0));
  # make sure the dimensions are in the right order
  $sum = $sum->xchg(0,1);
  # multiply with particle size distribution
  my $weighted = inner($sum/$xm/(re($ka)**2),$dQ2);
  # Divide by Frequency
  my $c_dispersion = 1/(1/$c-3*$phi*$weighted/$omega);
  my $dummy = (3*$phi*$weighted/$omega);
  #printf $dummy . "\n";
  return $c_dispersion;
}



sub alpha_Str_Faran {
  my ($f, $xm, $dQ2, $c, $h_L, $h_T, $g, $phi)=@_;
  # Calculate the wave numbers
  my $pi = 2*asin(1);
  # nächste Zeile- Wichtig für transition real -> complex
  my $ka=r2C(outer($xm,$f)*1e6/$c*$pi); # matmult ?
  my $ka_PL = $ka/$h_L;
  my $ka_PT = $ka/$h_T;
  #Calculate required length of Bessel function series
  my $n_b1 = n_B($ka);
  my $n_max = max($n_b1);
  my $n_b2 = n_B($ka_PL);
  my $n_max2 = max($n_b2);
  if ($n_max2 > $n_max) { $n_max = $n_max2;}
  my $n_b3 = n_B($ka_PT);
  my $n_max3 = max($n_b3);
  if ($n_max3 > $n_max) { $n_max = $n_max3;}
  # Calculate A_n (Kräuter (Diss), Hay & Mercer (1975) )
  my $An = A_n($ka ,$ka_PL ,$ka_PT ,$n_b1 ,$n_b2 ,$n_b3 ,$n_max ,$g);
  # A_n_real ?
  # Create Sum
  # 2*n+1
  my $n=sequence($n_max);
  $n = 2*$n+1;
  $n = $n->dummy(0,$ka->getdim(2));
  $n = $n->dummy(0,$ka->getdim(1));
  # (2n+1)*(abs(R_n))²
  my $Rn2 = (re($An)**2 + im($An)**2)*$n;
  my $re_Rn2=re($An);
  my $im_Rn2=im($An);
  #printf("re_Rn2 $re_Rn2\nim_Rn2 $im_Rn2\n");
  my $maxdim=$n_max-1;
  $Rn2 = $Rn2->slice(":,:,0:$maxdim");
  # Sum the above
  my $sum=dsumover($Rn2->xchg(2,0));
  # make sure the dimensions are in the right order
  $sum = $sum->xchg(0,1);
  # create alpha
  my $alpha = 3*$phi*$sum/(re($ka)**2);
  # multiply with particle size distribution
  my $weighted = inner($alpha/$xm,$dQ2);
  #printf("weighted $weighted\n");
  # Divide by Frequency
  my $alpha_div_f = $weighted/ $f;
  #printf("div f $alpha_div_f\n");
  # Result is in Np/m/MHz
  return $alpha_div_f;
}

sub alpha_Str_Faran_noweight {
  my ($f, $xm, $c, $h_L, $h_T, $g, $phi)=@_;
  # Calculate the wave numbers
  my $pi = 2*asin(1);
  my $ka=outer($xm,$f)*1e6/$c*$pi; # matmult ?
  my $ka_PL = $ka/$h_L;
  my $ka_PT = $ka/$h_T;
  #Calculate required length of Bessel function series
  my $n_b1 = n_B($ka);
  my $n_max = max($n_b1);
  my $n_b2 = n_B($ka_PL);
  my $n_max2 = max($n_b2);
  if ($n_max2 > $n_max) { $n_max = $n_max2;}
  my $n_b3 = n_B($ka_PT);
  my $n_max3 = max($n_b3);
  if ($n_max3 > $n_max) { $n_max = $n_max3;}
  # Calculate A_n (Kräuter (Diss), Hay & Mercer (1975) )
  my $An = A_n($ka ,$ka_PL ,$ka_PT ,$n_b1 ,$n_b2 ,$n_b3 ,$n_max ,$g);
  # Create Sum
  # 2*n+1
  my $n=sequence($n_max);
  $n = 2*$n+1;
  $n = $n->dummy(0,$ka->getdim(1));
  $n = $n->dummy(0,$ka->getdim(0));
  # (2n+1)*(abs(R_n))²
  my $Rn2 = (re($An)**2 + im($An)**2)*$n;
  my $re_Rn2=re($An);
  my $im_Rn2=im($An);
  #printf("re_Rn2 $re_Rn2\nim_Rn2 $im_Rn2\n");
  my $maxdim=$n_max-1;
  $Rn2 = $Rn2->slice(":,:,0:$maxdim");
  # Sum the above
  my $sum=dsumover($Rn2->xchg(2,0));
  # make sure the dimensions are in the right order
  $sum = $sum->xchg(0,1);
  # create alpha
  my $alpha = 3*$phi*$sum/($ka**2);
  # Divide by Frequency
  my $alpha_div_f = transpose(transpose($alpha/$xm)/$f);
  # printf("div f $alpha_div_f\n");
  # Result is in Np/m/MHz
  return $alpha_div_f;
}



# Calculates the scattering cross section for a
# combination of f and x
# for each single particle (no weighting)
sub A_Str_Faran {
  my ($f, $xm, $c, $h_L, $h_T, $g)=@_;
  # Calculate the wave numbers
  my $pi = 2*asin(1);
  my $ka=r2C(outer($xm,$f)*1e6/$c*$pi);
  my $ka_PL = $ka/$h_L;
  my $ka_PT = $ka/$h_T;
  my $k=2*$f*1e6/$c*$pi;
  #Calculate required length of Bessel function series
  my $n_b1 = n_B($ka);
  my $n_max = max($n_b1);
  my $n_b2 = n_B($ka_PL);
  my $n_max2 = max($n_b2);
  if ($n_max2 > $n_max) { $n_max = $n_max2;}
  my $n_b3 = n_B($ka_PT);
  my $n_max3 = max($n_b3);
  if ($n_max3 > $n_max) { $n_max = $n_max3;}
  # Calculate A_n (Kräuter (Diss), Hay & Mercer (1975) )
  my $An = A_n($ka, $ka_PL, $ka_PT, $n_b1, $n_b2, $n_b3, $n_max, $g);
  # Create Sum
  # 2*n+1
  my $n=sequence($n_max);
  $n = 2*$n+1;
  $n = $n->dummy(0,$ka->getdim(2));
  $n = $n->dummy(0,$ka->getdim(1));
  # (2n+1)*(abs(R_n))²
  my $Rn2 = (re($An)**2 + im($An)**2)*$n;
  my $maxdim=$n_max-1;
  $Rn2 = $Rn2->slice(":,:,0:$maxdim");
  #pdims($Rn2);
  # Sum the above
  my $sum=dsumover($Rn2->xchg(2,0));
  # make sure the dimensions are in the right order
  $sum = $sum->xchg(0,1);
  # create A
  my $A = 4*$pi*$sum/(transpose($k)**2);
  return $A;
}


