package PKM;
require 5.000;
require Exporter;
use Carp;

# PKM.pm
# (c) 2003 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.

# This module performs calculations of acoustic attenuation based
# upon the Coupled Phase Model (Dukhin, Goetz - Dispersion Technology).
# The code follows the calculations performed by a MathCad worksheet
# developed by Frank Babick (TU Dresden)
# Nomenclature and calculation was taken from the worksheet and
# adapted where needed to the specifics of Perl.

# The most important modification is the application of the PDL
# which allows for fast calculations of large arrays which are
# associated with particle size and frequency distributions

@ISA = qw(Exporter);
@EXPORT = qw(alpha_visc alpha_visc_noweight alpha_therm);

#!/usr/bin/perl

use strict;
use warnings;

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

my $dB=log(10)/20;

# Calculation of viscous attenuation based on given
# particle size distribution, frequency array
# and material properties

# Input (along with required units)
# f             : Frequencies vector (MHz)
# a             : Particle sizes vector (m)
# dQ3           : size distribution vector (in %, Sum(dQ3)= cumulative distribution)
# phi           : Solids volume fraction (0..1)
# K_P, K_F      : Modulus of Compressibility of _P_article and _F_luid (Pa)
# eta_F         : Dynamic _F_luid Viscosity (Pa*s)
# rho_F, rho_P  : Density of Fluid and Particle (kg/m³)

# Output (Units)
# attenuation alpha (dB/cm/Mhz) - vector

sub alpha_visc{
  my ($f,$a,$dQ3,$phi,$K_P, $K_F, $rho_P, $rho_F, $eta_F)=@_;
  my $C;
  my $denominator;
  my $KS;
  my $k_vis;
  my $alpha_vis;
  my $Sum1=0;
  my $Sum2=0;
  my $omega=omega($f);
  for(my $k=0;$k<$a->getdim(0);$k++) {
    $C=C_R($f,$phi,$a->slice("$k"),$rho_F, $eta_F);
    $denominator= ($C - $omega*$rho_P*i);
    $Sum1 = $Sum1+(($phi*$C*$dQ3->slice("$k")) / $denominator);
    $Sum2 = $Sum2+(($phi*$dQ3->slice("$k")) / $denominator);
  }
  $KS=effectiveK_S($phi, $K_P, $K_F );
  $k_vis=$omega/sqrt($KS) * sqrt( ((1-$phi)*$rho_F + $rho_P*$Sum1 ) /
     ((1-$phi+$Sum1)**2 -$omega*$Sum2*i*( (1-$phi)*$rho_F+$rho_P*$Sum1)));
  $alpha_vis=im($k_vis) / ( $f)/$dB / 100; # Units now: dB / MHz / cm
#  use PDL::Graphics2D;
#  use PDL::Graphics::PGPLOT;
#  dev('/XSERVE');
  #line($f,$alpha_vis, {LINESTYLE=>1}); # solid line
  return($alpha_vis);
}

# Calculation of viscous attenuation based on given
# particle size vector + frequency vector
# and material properties

# Input (along with required units)
# f             : Frequencies vector (MHz)
# a             : Particle sizes vector (m)
# phi           : Solids volume fraction (0..1)
# K_P, K_F      : Modulus of Compressibility of _P_article and _F_luid (Pa)
# eta_F         : Dynamic _F_luid Viscosity (Pa*s)
# rho_F, rho_P  : Density of Fluid and Particle (kg/m³)

# Output (Units)
# attenuation alpha (dB/cm/Mhz) - matrix (of order f,a)
# result is not weighted with a size distribution!

sub alpha_visc_noweight{
  my ($f,$a,$phi,$K_P, $K_F, $rho_P, $rho_F, $eta_F)=@_;
  my $C;
  my $denominator;
  my $KS;
  my $k_vis;
  my $alpha_vis;
  my $Sum1=zeroes($a->getdim(0));
  my $Sum2=zeroes($a->getdim(0));
  my $omega=omega($f);
  for(my $k=0;$k<$a->getdim(0);$k++) {
    $C=C_R($f,$phi,$a->slice("$k"),$rho_F, $eta_F);
    $denominator= ($C - $omega*$rho_P*i);
    $Sum1->slice("$k").=($phi*$C / $denominator);
    $Sum2->slice("$k").=($phi / $denominator);
  }
  #printf $Sum1 . "\n";
  $KS=effectiveK_S($phi, $K_P, $K_F );
  $k_vis=$omega/sqrt($KS) * sqrt( ((1-$phi)*$rho_F + $rho_P*transpose($Sum1) ) /
     ((1-$phi+transpose($Sum1))**2 -$omega*transpose($Sum2)*i*( (1-$phi)*$rho_F+$rho_P* transpose($Sum1))));
  $alpha_vis=im($k_vis) / ( $f)/$dB / 100; # Units now: dB / MHz / cm
  return(transpose($alpha_vis));
  # transpose: results in the matrix being of length(a) for 1st dimension,
  # length(f) for 2nd dimensions
  # required to conform with convention set by Faran implementation
}


# required functions for viscous attenuation


# effective module of compression
sub effectiveK_S {
  my ($phi, $K_P, $K_F)=@_;
  my $ret=($phi*($K_P**(-1)) + (1-$phi)*($K_F**(-1)) )**(-1);
  return $ret;
}

# effective density
sub effectiverho_S {
  my ($phi,$rho_P,$rho_F)=@_;
  my $ret=$rho_P*$phi+$rho_F*(1-$phi);
  return $ret;
}

# other functions
sub hzz { my ($x) =@_;
    my $ret=(-$x+i*(1+$x))/(($x)**2);
    return $ret;
}

sub omega { my ($f)= @_;
  my $pi = 2 *asin(1);
  my $ret=2*$pi*($f * 1e6); # MHz ->Hz

  return $ret;
}

sub I11 { my ($x) =@_;
    my $ret=i/-$x;
    return $ret;
}

sub I22 { my ($x) =@_;
    my $ret=-(3/2*(1+$x)+i*($x**2+3/2*$x));
    return $ret;
}


# friction coefficient
sub C_R {
  my ($f, $phi, $a, $rho_F, $eta_F)= @_;
  my $omega=omega($f);
  my $A=sqrt(($a**2*$omega*$rho_F)/2/$eta_F);
  my $b=$a*(1/$phi)**(1/3);
  my $kappa=$b/$a;
  my $B=$kappa*$A;
  my $h1=hzz($B)*exp((1-$kappa)*$A*(1+i));
  my $h2=hzz(-$B);
  my $test=hzz(-4);
  my $I1=I11($B)*exp((1-$kappa)*$A*(1+i))-I11($A);
  my $I2=-I11(-$B)+I11(-$A)*exp((1-$kappa)*$A*(1+i));
  my $I13=(I22($B)*exp((1-$kappa)*$A*(1+i))-I22($A))*($B**(-3));
  my $I23=-(-I22(-$B)+I22(-$A)*exp((1-$kappa)*$A*(1+i)))*($B**(-3));
  # Determinante
  my $Den=$I1*$h2*1 + $h1*$I23*1 + $I13*$I2*2*($B**-1) - (1*$h2*$I13 + 2*($B**-1)*$I23*$I1 + 1*$I2*$h1);
  my $C1=($h2-2*$I23/$B)*($Den**-1);
  my $C2=(-$h1+2*$I13/$B)*($Den**-1);
  my $H1=1/($A**2)*(3*i*(1/$A+1)-(3+2*$A));
  my $H2=exp((1-$kappa)*$A*(1+i))/($A**2)*(3*i*(1/$A-1)+(3-2*$A));
  my $Omega=($A**2)/3*($C1*$H1+$C2*$H2-4*i/3);
  my $res=Cconj($Omega)*9*$eta_F/2/($a**2);
#    printf("A = $A\n");
#    printf("kappa = $kappa\n");
#    printf("B = $B\n");
#    printf("h1 = $h1\n");
#    printf("h2 = $h2\n");
#    printf("I1 = $I1\n");
#    printf("I2 = $I2\n");
#    printf("I13 = $I13\n");
#    printf("I23 = $I23\n");
#    printf("Den = $Den\n");
#    printf("C1 = $C1\n");
#    printf("C2 = $C2\n");
#    printf("H1 = $H1\n");
#    printf("H2 = $H2\n");
#    printf("Omega = $Omega\n");
#    printf("res = $res\n");
#    printf("a $a\n");
#    printf("eta_f $eta_F\n");
  return $res;
}

# Calculation of the thermal attenuation

# Input (along with required units)
# f             : Frequencies vector (MHz)
# a             : Particle radius vector (m)
# dQ3           : size distribution vector (in %, Sum(dQ3)= cumulative distribution)
#                 needs to have the same length as vector a
# phi           : Solids volume fraction (0 ... 1)
# rho_F, rho_P  : Density of _F_luid and _P_article (kg/m³)
# K_F           : Modulus of Compressibility of the _F_luid (Pa)
# c_PF, c_PP    : specific heat capacity (W*s/kg/K)
# tau_F, tau_P  : Heat transfer coefficient (W/K/m)
# eta_F         : Dynamic _F_luid Viscosity (Pa*s)
# beta_P, beta_F: isobaric expansion coefficient of Particle and Fluid (1/K)

# Output (Units)
# attenuation alpha (dB/cm/Mhz) - vector

sub alpha_therm {
  my ($f,$a,$dQ3,$phi, $rho_P, $rho_F, $c_PF, $c_PP, $tau_F, $tau_P, $eta_F, $K_F, $beta_P,$beta_F)=@_;
  # reference temperature (K)
  my $T_0 = 298;
  # circular frequency
  my $omega_F = omega($f);
  # Velocity of sound (compressional wave, m/s)
  my $c_KF = sqrt($K_F/$rho_F);
  # adiabatic exponent of the fluid
  my $gamma_F=1+($beta_F**2 *$K_F* $T_0)/($c_PF*$rho_F);
  # TF is a helper variable
  my $TF=($gamma_F -1) *
         (1-($beta_P*$rho_F*$c_PF)/($beta_F*$rho_P*$c_PP))**2;
  # dimensionless wave numbers (Fluid and Particle)
  my $z_thF =z_th($f,$a,$rho_F,$c_PF,$tau_F);
  my $z_thP =z_th($f,$a,$rho_P,$c_PP,$tau_P);
  # Calculate tangent explicitly - speeds up calculation of next step
  # Ctan is used because z_thP is complex
  my $Ctan_z = Ctan($z_thP);
  # H is another helper variable
  my $H = 1/( 1 /(1-(i)*$z_thF)
            -$tau_F/$tau_P
             * $Ctan_z/($Ctan_z-$z_thP));
  # thermal parameter
  my $X_F = X($f,$a,$rho_F,$c_PF,$tau_F);
  # combine for each particle size and weight according to size distribution
  my $term1 = re($H)/($X_F**2) * $dQ3;
  # sum over all particle sizes, gives a vector of the same size as the
  # frequency vector
  my $sum = sumover($term1);
  my $res = 3*$phi/4/$c_KF*$TF*  $omega_F * $sum;
  # some debugging: output intermediate calculation results
#  printf "re(H)" . re($H)->slice("0,:") . "\n";
#  printf "X_F" . $X_F->slice("0,:") . "\n";
#  printf "omega" . $omega_F . "\n";
#  printf "re(H)/(X_F**2) * dQ3" . $term1->slice(":,1") . "\n";
#  printf "sumover(re(H)/(X_F**2)*dQ3) ". $sum . "\n";
  my $res_dB_cm_Mhz = $res/ ( $f)/$dB / 100;
  # return result
  # units are db/MHz/cm
  return $res_dB_cm_Mhz;
}

# required functions for alpha_therm
# thermal parameter

sub X {
  my ($f,$a,$rho,$c_p,$tau) = @_;
  #my $res=sqrt(matmult($a**2,omega($f))* $rho*$c /2 / $tau);
  my $res=sqrt(outer($a**2,omega($f))* $rho*$c_p /2 / $tau);
  return $res;
}

# dim.-less thermal wave number
sub z_th {
  my ($f,$a,$rho,$c,$tau) = @_;
  my $X_=X($f,$a,$rho,$c,$tau);
  my $res=(1+i)*$X_;
  return $res;
}


1;
