#!/usr/bin/env perl
#;-*- Perl -*-

# Created by dano 6-19-09
# LAST MODIFIED by dano: 6-4-08
# Reads cube files (electrostatic potential) returns value at r_ions
######################################################
# Still need to add periodic boundy cond
######################################################

use FindBin qw($Bin);
use lib "$Bin";
use Vasp;
use Math::Trig;
$fact = 180/pi;

print "\n";
@args = @ARGV;
(@args >= 1 || @args >= 2) || die "usage: cubevalue.pl <cube file> <method>\n";
$inputfilename = $args[0];
if(@args == 2){
    $method = @args[1];
}else{
    $method = 1;
}

# Assume atomic units 
$unittype = "bohr";
$ang2bohr_units = 1.889725992; # ang--->bohr
$eV2hartree_units = 0.036749; # eV---> hartree
$bohr2ang_units = .529177248; # bohr--->ang
$hartree2eV_units = 27.2116248; # hartree ---> eV
$unit_flag = 1;
print "     USING '$unittype' FOR OUTPUT UNITS\n";

# prompt USER for file format
#print "WAS CUBE FILE GENERATED BY chg2cube.pl? 1=YES or 2=NO:  ";
#chomp($filetype_ = <STDIN>);
$filetype_ = 2;
if ($filetype_ == 1) {$filetype= "vasp";}
if ($filetype_ == 2) {$filetype= "cube";}
if ($filetype_ == 1 or $filetype_ == 2) {
#  print "     READING '$filetype' FILE FORMAT\n";
    print "\n";
} else {
    die "     YOU MUST ENTER 1 or 2 FOR FILE CONVERSION INFO\n";
}

# READ FILE
$inputfile = "";
open (IN,$inputfilename);
while (<IN>) {
    $_ =~ s/^\s+//g;
    $inputfile .= $_;
}
close (IN);

@inputfile=split(/\n/,$inputfile);

#skip first two lines 
$line_ = $inputfile[2];
$line_ =~ s/^\s+//;
@line = split(/\s+/,$line_);
$total_atoms = $line[0];
$origin->[0] = $line[1];
$origin->[1] = $line[2];
$origin->[2] = $line[3];

##### read in grid and lattice ##########
for ($i=0; $i<3; $i++) {
    $line = $inputfile[$i+3];
    $line =~ s/^\s+//;
    @line = split(/\s+/,$line);
    # first value is # of voxles in that direction
    $grid->[$i] = $line[0];
    # This is how cube reads in the basis
    for ($j=0; $j<3; $j++) {
        $basis->[$j][$i] = $line[$j+1];
    }
}

# find inverse basis
$inverse = &inverse($basis); # call function inverse from Vasp
#print "$inverse->[0][0] $inverse->[0][2] $inverse->[0][1]\n";

# calculate the volume of a cube (expansion by minors)
$dvol = $basis->[0][0]*($basis->[1][1]*$basis->[2][2] - $basis->[2][1]*$basis->[1][2])
      - $basis->[1][0]*($basis->[0][1]*$basis->[2][2] - $basis->[2][1]*$basis->[0][2])
      + $basis->[2][0]*($basis->[0][1]*$basis->[1][2] - $basis->[1][1]*$basis->[0][3]);
$dvol = abs($dvol);
$grid_points = $grid->[0]*$grid->[1]*$grid->[2];
$box_vol = $dvol*$grid_points;

# scale basis to size of box!
for ($i=0; $i<3; $i++) {
    for ($j=0; $j<3; $j++) {
        $scale_basis->[$j][$i] = $grid->[$i]*$basis->[$j][$i];
    }
}

##### read in the atoms positions (are in cart coordinates)#####
$index = 6;
for ($i=0; $i<$total_atoms; $i++) {
    $line = $inputfile[$i+$index];
    $line =~ s/^\s+//;
    @line = split(/\s+/,$line);  
    $atom_type->[$i] = $line[0];
    for ($j=2; $j<5; $j++) {
        $coordinates->[$i][$j-2] = $line[$j];
    }
}

##### read in grid data ##########
$index = $index+$total_atoms;
$data_lines = int($grid_points/6+1);
$ii = 0;
for ($i=0; $i<$data_lines; $i++) {
#for ($i=0; $i<3; $i++) {
    $line = $inputfile[$i+$index];
    $line =~ s/^\s+//;
    @line = split(/\s+/,$line);
    for ($j=0; $j<@line; $j++) {
    # for cube files density is z-inner, y-middle, x-outter loop
        $density->[$ii] = $line[$j];
        $ii++;
    }
}

# it would be nice if grid data was in t x,y,z format
$ii = 0;
for ($i=0; $i<$grid->[0]; $i++) {
    for ($j=0; $j<$grid->[1]; $j++) {
        for ($k=0; $k<$grid->[2]; $k++) {
            $rho->[$i][$j][$k] = $density->[$ii];
            $ii++;
        }
    }
}

# find value of cube-file at posions of ions
for ($i=0; $i<$total_atoms; $i++) {
    # change xyz to grid value (function)
    $ijk = &xyz2grid($coordinates->[$i][0],$coordinates->[$i][1],$coordinates->[$i][2]);
    #  print "$ijk->[0] $ijk->[1] $ijk->[2]\n";
    #interpolate for density value
    if ($method == 1){
        $rho_R->[$i] = &lineinterp($ijk->[0],$ijk->[1],$ijk->[2]);}
    if ($method == 2){
        $rho_R->[$i] = &quadinterp($ijk->[0],$ijk->[1],$ijk->[2]);}
    if ($method == 3){
        $rho_R->[$i] = &cubicinterp($ijk->[0],$ijk->[1],$ijk->[2]);}
    if ($method == 4){
        $rho_R->[$i] = &minvalue($ijk->[0],$ijk->[1],$ijk->[2]);}
    #  print "$rho_R->[$i]\n";
}

########## output values at atom posions ########################
if(@args >= 2) {
    $outputfilename = $args[1];
} else {
    $outputfilename = "out.chempot";
}
open (OUT,">$outputfilename");

print OUT "Cartesian Coordinates (Bohr)\n";
print OUT "   Zval       x            y            z        u(x,y,z) hartree\n";
# print data
print "\n";
print "Volume DATA R_ions (hartree)\n";
for ($i=0;$i<$total_atoms;$i++) {
    printf OUT "%5i %12.6f %12.6f %12.6f %12.6f\n",$atom_type->[$i],$coordinates->[$i][0],$coordinates->[$i][1],$coordinates->[$i][2],$rho_R->[$i];
    #  print "Volume DATA R_ions (hartree)\n"
    print "  $rho_R->[$i]\n";
}

print OUT "\n";
print OUT "\n";
print OUT "Cartesian Coordinates (Angstroms)\n";
print OUT "   Zval       x            y            z        u(x,y,z) eV\n";
# print data
print "\n";
print "Volume DATA R_ions (eV)\n";
for ($i=0; $i<$total_atoms; $i++) {
    $rho_R_eV->[$i] = $rho_R->[$i]*$hartree2eV_units;
    printf OUT "%5i %12.6f %12.6f %12.6f %12.6f\n",$atom_type->[$i],$coordinates->[$i][0]*$bohr2ang_units,$coordinates->[$i][1]*$bohr2ang_units,$coordinates->[$i][2]*$bohr2ang_units,$rho_R_eV->[$i];
    print "  $rho_R_eV->[$i]\n";
}
print OUT "\n";
close(OUT);
print "DONE!\n";

#############################################################################
# subroutine to find xyz coordinates of point on the grid
#############################################################################
# assign cartensin coordiantes to grid points (call it center of the volume)
sub grid2xyz {
    my($i,$j,$k) = @_;
    my($xyz);
    if ($filetype_ == 1) {
        $xyz->[0] = (0.5+$i)*$basis->[0][0] + (0.5+$j)*$basis->[0][1] + (0.5+$k)*$basis->[0][2];
        $xyz->[1] = (0.5+$i)*$basis->[1][0] + (0.5+$j)*$basis->[1][1] + (0.5+$k)*$basis->[1][2];
        $xyz->[2] = (0.5+$i)*$basis->[2][0] + (0.5+$j)*$basis->[2][1] + (0.5+$k)*$basis->[2][2];
    }
    if ($filetype_ == 2) {
        $xyz->[0] = ($i)*$basis->[0][0] + ($j)*$basis->[0][1] + ($k)*$basis->[0][2];
        $xyz->[1] = ($i)*$basis->[1][0] + ($j)*$basis->[1][1] + ($k)*$basis->[1][2];
        $xyz->[2] = ($i)*$basis->[2][0] + ($j)*$basis->[2][1] + ($k)*$basis->[2][2];
    }
    return($xyz);
}

#############################################################################
# subroutine to exact grid point for xyz input (does not return int)
#############################################################################
# (call it center of the volume of grid)
sub xyz2grid {
    my($x,$y,$z) = @_;
    my($ijk);
    $ijk->[0] = $x*$inverse->[0][0] + $y*$inverse->[1][0] + $z*$inverse->[2][0];
    $ijk->[1] = $x*$inverse->[0][1] + $y*$inverse->[1][1] + $z*$inverse->[2][1];
    $ijk->[2] = $x*$inverse->[0][2] + $y*$inverse->[1][2] + $z*$inverse->[2][2];
    if ($filetype_ == 1) {
        # generated from vasp format
        $ijk->[0] -= 0.5;
        $ijk->[1] -= 0.5;
        $ijk->[2] -= 0.5;
    }
    return($ijk);
}

#############################################################################
# linear interpolation for finding density value on grid
#############################################################################
sub lineinterp {
    my($i,$j,$k) = @_;
    # gridpoint prev and next
    my($ip, $in, $jp, $jn, $kp, $kn);
    # distance from prev gridpoint
    my($di, $dj, $dk);
    # interpolation matrix
    my($interp, $val);

    # check to see if outside grid bounries
    if ($i>0 and $i+1<$grid->[0]) {
        $ip = int($i);
        $in = int($i+1);
        $di = $i-$ip;
    } else {
        $ip = $grid->[0];
        $in=0;
        if ($i>0) {
            $di = 1+$i;
        } else {
            $di = $i-$ip;
        }
    }
    if ($j>0 and $j+1<$grid->[1]) {
        $jp = int($j);
        $jn = int($j+1);
        $dj = $j-$jp;
    } else {
        $jp = $grid->[1];
        $jn = 0;
        if ($j>0) {
            $dj = 1+$j;
        } else {
            $dj = $j-$jp;
        }
    }
    if ($k>0 and $k+1<$grid->[2]) {
        $kp = int($k);
        $kn = int($k+1);
        $dk = $k-$kp;
    } else {
        $kp = $grid->[2];
        $kn = 0;
        if ($k>0) {
            $dk = 1+$k;
        } else {
            $dk = $k-$kp;
        }
    }

    # get 8 values around the real
    $interp->[0][0][0] = $rho->[$ip][$jp][$kp];
    $interp->[2][0][0] = $rho->[$in][$jp][$kp];
    $interp->[0][2][0] = $rho->[$ip][$jn][$kp];
    $interp->[2][2][0] = $rho->[$in][$jn][$kp];
    $interp->[0][0][2] = $rho->[$ip][$jp][$kn];
    $interp->[2][0][2] = $rho->[$in][$jp][$kn];
    $interp->[0][2][2] = $rho->[$ip][$jn][$kn];
    $interp->[2][2][2] = $rho->[$in][$jn][$kn];

    ##### interpolate 12 edges #######
    #$interp->[1][0][0] = (1-$di)*$interp->[0][0][0] + $di*$interp->[2][0][0]; 
    $interp->[0][1][0] = (1-$dj)*$interp->[0][0][0] + $dj*$interp->[0][2][0]; 
    #$interp->[0][0][1] = (1-$dk)*$interp->[0][0][0] + $dk*$interp->[0][0][2];

    $interp->[2][1][0] = (1-$dj)*$interp->[2][0][0] + $dj*$interp->[2][2][0]; 
    #$interp->[2][0][1] = (1-$dk)*$interp->[2][0][0] + $dk*$interp->[2][0][2]; 

    #$interp->[1][2][0] = (1-$di)*$interp->[0][2][0] + $di*$interp->[2][2][0]; 
    #$interp->[0][2][1] = (1-$dk)*$interp->[0][2][0] + $dk*$interp->[0][2][2]; 

    #$interp->[1][0][2] = (1-$di)*$interp->[0][0][2] + $di*$interp->[2][0][2]; 
    $interp->[0][1][2] = (1-$dj)*$interp->[0][0][2] + $dj*$interp->[0][2][2]; 

    #$interp->[1][2][2] = (1-$di)*$interp->[0][2][2] + $di*$interp->[2][2][2]; 
    $interp->[2][1][2] = (1-$dj)*$interp->[2][0][2] + $dj*$interp->[2][2][2]; 
    #$interp->[2][2][1] = (1-$dk)*$interp->[2][2][0] + $dk*$interp->[2][2][2];

    #### interpolate 6 faces ########
    #$interp->[1][1][0] = (1-$dj)*$interp->[1][0][0] + $dj*$interp->[1][2][0];
    #$interp->[1][0][1] = (1-$dk)*$interp->[1][0][0] + $dk*$interp->[1][0][2];
    $interp->[0][1][1] = (1-$dk)*$interp->[0][1][0] + $dk*$interp->[0][1][2];

    #$interp->[1][1][2] = (1-$dj)*$interp->[1][0][2] + $dj*$interp->[1][2][2];
    #$interp->[1][2][1] = (1-$dk)*$interp->[1][2][0] + $dk*$interp->[1][2][2];
    $interp->[2][1][1] = (1-$dk)*$interp->[2][1][0] + $dk*$interp->[2][1][2];

    #### interpolate to the center #####
    $interp->[1][1][1] = (1-$di)*$interp->[0][1][1] + $di*$interp->[2][1][1];
    $val = $interp->[1][1][1];
    return($val);
}

#############################################################################
# use the min value surrounding gridpoints 
#############################################################################
sub minvalue {
    my($i,$j,$k) = @_;
    # gridpoint prev and next
    my($ip, $in, $jp, $jn, $kp, $kn);
    # distance from prev gridpoint
    my($di, $dj, $dk);
    # interpolation matrix
    my($interp, $val);

    # check to see if outside grid boundaries
    if ($i>0 and $i+1<$fft_grid->[0]) {
        $ip = int($i);
        $in = int($i+1);
        $di = $i-$ip;
    } else {
        $ip = $fft_grid->[0];
        $in = 0;
        if ($i>0) {
            $di = 1+$i;
        } else {
            $di = $i-$ip;
        }
    }
    if ($j>0 and $j+1<$fft_grid->[1]) {
        $jp = int($j);
        $jn = int($j+1);
        $dj = $j-$jp;
    } else {
        $jp = $fft_grid->[1];
        $jn = 0;
        if ($j>0) {
            $dj = 1+$j;
        } else {
            $dj = $j-$jp;
        }
    }
    if ($k>0 and $k+1<$fft_grid->[2]) {
        $kp = int($k);
        $kn = int($k+1);
        $dk = $k-$kp;
    } else {
        $kp = $fft_grid->[2];
        $kn = 0;
        if ($k>0) {
            $dk = 1+$k;
        } else {
            $dk = $k-$kp;
        }
    }

    # get 8 values around the real
    $interp->[0][0][0] = $rho->[$ip][$jp][$kp];
    $interp->[2][0][0] = $rho->[$in][$jp][$kp];
    $interp->[0][2][0] = $rho->[$ip][$jn][$kp];
    $interp->[2][2][0] = $rho->[$in][$jn][$kp];
    $interp->[0][0][2] = $rho->[$ip][$jp][$kn];
    $interp->[2][0][2] = $rho->[$in][$jp][$kn];
    $interp->[0][2][2] = $rho->[$ip][$jn][$kn];
    $interp->[2][2][2] = $rho->[$in][$jn][$kn];

    $val = $interp->[0][0][0];
    if ($val>$interp->[2][0][0]) { $val = $interp->[2][0][0]; }
    if ($val>$interp->[0][2][0]) { $val = $interp->[0][2][0]; }
    if ($val>$interp->[2][2][0]) { $val = $interp->[2][2][0]; }
    if ($val>$interp->[0][0][2]) { $val = $interp->[0][0][2]; }
    if ($val>$interp->[2][0][2]) { $val = $interp->[2][0][2]; }
    if ($val>$interp->[0][2][2]) { $val = $interp->[0][2][2]; }
    if ($val>$interp->[2][2][2]) { $val = $interp->[2][2][2]; }

    return($val);
}
##############################################################################
#############################################################################
# quadradic interpolation for finding density value on grid
#############################################################################
sub quadinterp {
    my($i,$j,$k) = @_;
    # gridpoint prev and next
    my($ip, $jp, $kp);
    # distance from prev gridpoint
    my($di, $dj, $dk);
    # interpolation matrix
    my($interp);
    my($val);

    # find gridpoints around point
    if ($i-int($i) >= 0.5) {
        $ip = int($i);
    } else {
        $ip = int($i-1);
    }
    my $in = int($ip+1);
    my $i2n = int($ip+2);
    $di = $i-$ip;
    if ($j-int($j) >= 0.5) {
        $jp = int($j);
    } else {
        $jp = int($j-1);
    }
    my $jn = int($jp+1);
    my $j2n = int($jp+2);
    $dj = $j-$jp;
    if ($k-int($k) >= 0.5) {
        $kp = int($k);
    } else {
        $kp = int($k-1);
    }
    my $kn = int($kp+1);
    my $k2n = int($kp+2);
    $dk = $k-$kp;

    # get 27 values around the point of interest
    $interp->[0][0][0] = $rho->[$ip][$jp][$kp];
    $interp->[2][0][0] = $rho->[$in][$jp][$kp];
    $interp->[3][0][0] = $rho->[$i2n][$jp][$kp];
    $interp->[0][2][0] = $rho->[$ip][$jn][$kp];
    $interp->[2][2][0] = $rho->[$in][$jn][$kp];
    $interp->[3][2][0] = $rho->[$i2n][$jn][$kp];
    $interp->[0][3][0] = $rho->[$ip][$j2n][$kp];
    $interp->[2][3][0] = $rho->[$in][$j2n][$kp];
    $interp->[3][3][0] = $rho->[$i2n][$j2n][$kp];
    $interp->[0][0][2] = $rho->[$ip][$jp][$kn];
    $interp->[2][0][2] = $rho->[$in][$jp][$kn];
    $interp->[3][0][2] = $rho->[$i2n][$jp][$kn];
    $interp->[0][2][2] = $rho->[$ip][$jn][$kn];
    $interp->[2][2][2] = $rho->[$in][$jn][$kn];
    $interp->[3][2][2] = $rho->[$i2n][$jn][$kn];
    $interp->[0][3][2] = $rho->[$ip][$j2n][$kn];
    $interp->[2][3][2] = $rho->[$in][$j2n][$kn];
    $interp->[3][3][2] = $rho->[$i2n][$j2n][$kn];
    $interp->[0][0][3] = $rho->[$ip][$jp][$k2n];
    $interp->[2][0][3] = $rho->[$in][$jp][$k2n];
    $interp->[3][0][3] = $rho->[$i2n][$jp][$k2n];
    $interp->[0][2][3] = $rho->[$ip][$jn][$k2n];
    $interp->[2][2][3] = $rho->[$in][$jn][$k2n];
    $interp->[3][2][3] = $rho->[$i2n][$jn][$k2n];
    $interp->[0][3][3] = $rho->[$ip][$j2n][$k2n];
    $interp->[2][3][3] = $rho->[$in][$j2n][$k2n];
    $interp->[3][3][3] = $rho->[$i2n][$j2n][$k2n];

    ###### subroutine to fit a parabla to and interpolate y at x #####
    sub parabfit {
        #### intput is three x,f(x) points and value x for to solve for
        my($x1,$y1,$x2,$y2,$x3,$y3,$x) = @_;

        my $dx12 = ($x1 - $x2);
        my $sq_dx12 = ($x1**2 - $x2**2);
        my $dx13 = ($x1-$x3);
        my $sq_dx13 = ($x1**2 - $x3**2);
        my $dy12 = ($y1 - $y2);
        my $dy13 = ($y1 - $y3);

        my $b = ($dy13/$sq_dx13 - $dy12/$sq_dx12)/($dx13/$sq_dx13 - $dx12/$sq_dx12);
        my $a = ($dy12-$b*$dx12)/$sq_dx12;
        my $c = $y1 - $a*$x1**2 - $b*$x1;
        my $y = $a*$x**2 + $b*$x+$c;
        return($y);
    };

    ##### quadratically interpolate along the x direction #######
    $interp->[1][0][0] = &parabfit(0.,$interp->[0][0][0],1.,$interp->[2][0][0],2.,$interp->[3][0][0],$di);
    $interp->[1][2][0] = &parabfit(0.,$interp->[0][2][0],1.,$interp->[2][2][0],2.,$interp->[3][2][0],$di);
    $interp->[1][3][0] = &parabfit(0.,$interp->[0][3][0],1.,$interp->[2][3][0],2.,$interp->[3][3][0],$di);
    $interp->[1][0][2] = &parabfit(0.,$interp->[0][0][2],1.,$interp->[2][0][2],2.,$interp->[3][0][2],$di);
    $interp->[1][2][2] = &parabfit(0.,$interp->[0][2][2],1.,$interp->[2][2][2],2.,$interp->[3][2][2],$di);
    $interp->[1][3][2] = &parabfit(0.,$interp->[0][3][2],1.,$interp->[2][3][2],2.,$interp->[3][3][2],$di);
    $interp->[1][0][3] = &parabfit(0.,$interp->[0][0][3],1.,$interp->[2][0][3],2.,$interp->[3][0][3],$di);
    $interp->[1][2][3] = &parabfit(0.,$interp->[0][2][3],1.,$interp->[2][2][3],2.,$interp->[3][2][3],$di);
    $interp->[1][3][3] = &parabfit(0.,$interp->[0][3][3],1.,$interp->[2][3][3],2.,$interp->[3][3][3],$di);

    ##### quadratically interpolate along the y direction #######
    $interp->[1][1][0] = &parabfit(0.,$interp->[1][0][0],1.,$interp->[1][2][0],2.,$interp->[1][3][0],$dj);
    $interp->[1][1][2] = &parabfit(0.,$interp->[1][0][2],1.,$interp->[1][2][2],2.,$interp->[1][3][2],$dj);
    $interp->[1][1][3] = &parabfit(0.,$interp->[1][0][3],1.,$interp->[1][2][3],2.,$interp->[1][3][3],$dj);

    ##### quadratically interpolate along the z direction #######
    $interp->[1][1][1] = &parabfit(0.,$interp->[1][1][0],1.,$interp->[1][1][2],2.,$interp->[1][1][3],$dk);

    $val = $interp->[1][1][1];
    return($val);
}

#############################################################################
# cubic interpolation for finding density value on grid
#############################################################################
sub cubicinterp {
    my($i,$j,$k) = @_;
    # find gridpoints around point
    my $ip = int($i);
    my $i2p = int($ip-1);
    my $in = int($ip+1);
    my $i2n = int($ip+2);
    my $di = $i - $i2p;
    my $jp = int($j);
    my $j2p = int($jp-1);
    my $jn = int($jp+1);
    my $j2n = int($jp+2);
    # distance from prev gridpoint
    my $dj = $j - $j2p;
    my $kp = int($k);
    my $k2p = int($kp-1);
    my $kn = int($kp+1);
    my $k2n = int($kp+2);
    my $dk = $k - $k2p;
    # interpolation matrix
    my($interp);
    my($val);

  # get 48 values around the point of interest
    $interp->[0][0][0] = $rho->[$i2p][$j2p][$k2p];
    $interp->[1][0][0] = $rho->[$ip][$j2p][$k2p];
    $interp->[3][0][0] = $rho->[$in][$j2p][$k2p];
    $interp->[4][0][0] = $rho->[$i2n][$j2p][$k2p];
    $interp->[0][1][0] = $rho->[$i2p][$jp][$k2p];
    $interp->[1][1][0] = $rho->[$ip][$jp][$k2p];
    $interp->[3][1][0] = $rho->[$in][$jp][$k2p];
    $interp->[4][1][0] = $rho->[$i2n][$jp][$k2p];
    $interp->[0][3][0] = $rho->[$i2p][$jn][$k2p];
    $interp->[1][3][0] = $rho->[$ip][$jn][$k2p];
    $interp->[3][3][0] = $rho->[$in][$jn][$k2p];
    $interp->[4][3][0] = $rho->[$i2n][$jn][$k2p];
    $interp->[0][4][0] = $rho->[$i2p][$j2n][$k2p];
    $interp->[1][4][0] = $rho->[$ip][$j2n][$k2p];
    $interp->[3][4][0] = $rho->[$in][$j2n][$k2p];
    $interp->[4][4][0] = $rho->[$i2n][$j2n][$k2p];
    $interp->[0][0][1] = $rho->[$i2p][$j2p][$kp];
    $interp->[1][0][1] = $rho->[$ip][$j2p][$kp];
    $interp->[3][0][1] = $rho->[$in][$j2p][$kp];
    $interp->[4][0][1] = $rho->[$i2n][$j2p][$kp];
    $interp->[0][1][1] = $rho->[$i2p][$jp][$kp];
    $interp->[1][1][1] = $rho->[$ip][$jp][$kp];
    $interp->[3][1][1] = $rho->[$in][$jp][$kp];
    $interp->[4][1][1] = $rho->[$i2n][$jp][$kp];
    $interp->[0][3][1] = $rho->[$i2p][$jn][$kp];
    $interp->[1][3][1] = $rho->[$ip][$jn][$kp];
    $interp->[3][3][1] = $rho->[$in][$jn][$kp];
    $interp->[4][3][1] = $rho->[$i2n][$jn][$kp];
    $interp->[0][4][1] = $rho->[$i2p][$j2n][$kp];
    $interp->[1][4][1] = $rho->[$ip][$j2n][$kp];
    $interp->[3][4][1] = $rho->[$in][$j2n][$kp];
    $interp->[4][4][1] = $rho->[$i2n][$j2n][$kp];
    $interp->[0][0][3] = $rho->[$i2p][$j2p][$kn];
    $interp->[1][0][3] = $rho->[$ip][$j2p][$kn];
    $interp->[3][0][3] = $rho->[$in][$j2p][$kn];
    $interp->[4][0][3] = $rho->[$i2n][$j2p][$kn];
    $interp->[0][1][3] = $rho->[$i2p][$jp][$kn];
    $interp->[1][1][3] = $rho->[$ip][$jp][$kn];
    $interp->[3][1][3] = $rho->[$in][$jp][$kn];
    $interp->[4][1][3] = $rho->[$i2n][$jp][$kn];
    $interp->[0][3][3] = $rho->[$i2p][$jn][$kn];
    $interp->[1][3][3] = $rho->[$ip][$jn][$kn];
    $interp->[3][3][3] = $rho->[$in][$jn][$kn];
    $interp->[4][3][3] = $rho->[$i2n][$jn][$kn];
    $interp->[0][4][3] = $rho->[$i2p][$j2n][$kn];
    $interp->[1][4][3] = $rho->[$ip][$j2n][$kn];
    $interp->[3][4][3] = $rho->[$in][$j2n][$kn];
    $interp->[4][4][3] = $rho->[$i2n][$j2n][$kn];
    $interp->[0][0][4] = $rho->[$i2p][$j2p][$k2n];
    $interp->[1][0][4] = $rho->[$ip][$j2p][$k2n];
    $interp->[3][0][4] = $rho->[$in][$j2p][$k2n];
    $interp->[4][0][4] = $rho->[$i2n][$j2p][$k2n];
    $interp->[0][1][4] = $rho->[$i2p][$jp][$k2n];
    $interp->[1][1][4] = $rho->[$ip][$jp][$k2n];
    $interp->[3][1][4] = $rho->[$in][$jp][$k2n];
    $interp->[4][1][4] = $rho->[$i2n][$jp][$k2n];
    $interp->[0][3][4] = $rho->[$i2p][$jn][$k2n];
    $interp->[1][3][4] = $rho->[$ip][$jn][$k2n];
    $interp->[3][3][4] = $rho->[$in][$jn][$k2n];
    $interp->[4][3][4] = $rho->[$i2n][$jn][$k2n];
    $interp->[0][4][4] = $rho->[$i2p][$j2n][$k2n];
    $interp->[1][4][4] = $rho->[$ip][$j2n][$k2n];
    $interp->[3][4][4] = $rho->[$in][$j2n][$k2n];
    $interp->[4][4][4] = $rho->[$i2n][$j2n][$k2n];

    ###### subroutine to fit a cubic to 4 points and interpolate y at x #####
    sub cubefit {
        #### intput is three x,f(x) points and value x for to solve for
        my($x1,$y1,$x2,$y2,$x3,$y3,$x4,$y4,$x) = @_;

        my $dx12 = ($x1-$x2);
        my $sq_dx12 = ($x1**2-$x2**2);
        my $cu_dx12 = ($x1**3-$x2**3);
        my $dx13 = ($x1-$x3);
        my $sq_dx13 = ($x1**2-$x3**2);
        my $cu_dx13 = ($x1**3-$x3**3);
        my $dx14 = ($x1-$x4);
        my $sq_dx14 = ($x1**2-$x4**2);
        my $cu_dx14 = ($x1**3-$x4**3);
        my $dy12 = ($y1-$y2);
        my $dy13 = ($y1-$y3);
        my $dy14 = ($y1-$y4);

        my $a = ($dy12*$dx13 - $dy13*$dx12)*($sq_dx12*$dx14 - $sq_dx14*$dx12);
        $a = $a - ($dy12*$dx14 - $dy14*$dx12)*($sq_dx12*$dx13 - $sq_dx13*$dx12);
        my $a_dnom = ($cu_dx12*$dx14 - $cu_dx14*$dx12)*($sq_dx12*$dx13 - $sq_dx13*$dx12);
        $a_dnom = $a_dnom - ($cu_dx12*$dx13 - $cu_dx13*$dx12)*($sq_dx12*$dx14 - $sq_dx14*$dx12);
        $a = $a/$a_dnom;
        my $b = $dy12*$dx13 - $dy13*$dx12 - $a*($cu_dx12*$dx13 - $cu_dx13*$dx12);
        $b = $b/($sq_dx12*$dx13 - $sq_dx13*$dx12);
        my $c = ($dy12 - $a*$cu_dx12 - $b*$sq_dx12)/$dx12;
        my $d = $y1 - $a*$x1**3 - $b*$x1**2 - $c*$x1;
        my $y = $a*$x**3 + $b*$x**2 + $c*$x + $d;
        return($y);
    };

    ##### cubically interpolate along the x direction #######
    $interp->[2][0][0] = &cubefit(0.,$interp->[0][0][0],1.,$interp->[1][0][0],2.,$interp->[3][0][0],3.,$interp->[3][0][0],$di);
    $interp->[2][1][0] = &cubefit(0.,$interp->[0][1][0],1.,$interp->[1][1][0],2.,$interp->[3][1][0],3.,$interp->[3][1][0],$di);
    $interp->[2][3][0] = &cubefit(0.,$interp->[0][3][0],1.,$interp->[1][3][0],2.,$interp->[3][3][0],3.,$interp->[3][3][0],$di);
    $interp->[2][4][0] = &cubefit(0.,$interp->[0][4][0],1.,$interp->[1][4][0],2.,$interp->[3][4][0],3.,$interp->[3][4][0],$di);
    $interp->[2][0][1] = &cubefit(0.,$interp->[0][0][1],1.,$interp->[1][0][1],2.,$interp->[3][0][1],3.,$interp->[3][0][1],$di);
    $interp->[2][1][1] = &cubefit(0.,$interp->[0][1][1],1.,$interp->[1][1][1],2.,$interp->[3][1][1],3.,$interp->[3][1][1],$di);
    $interp->[2][3][1] = &cubefit(0.,$interp->[0][3][1],1.,$interp->[1][3][1],2.,$interp->[3][3][1],3.,$interp->[3][3][1],$di);
    $interp->[2][4][1] = &cubefit(0.,$interp->[0][4][1],1.,$interp->[1][4][1],2.,$interp->[3][4][1],3.,$interp->[3][4][1],$di);
    $interp->[2][0][3] = &cubefit(0.,$interp->[0][0][3],1.,$interp->[1][0][3],2.,$interp->[3][0][3],3.,$interp->[3][0][3],$di);
    $interp->[2][1][3] = &cubefit(0.,$interp->[0][1][3],1.,$interp->[1][1][3],2.,$interp->[3][1][3],3.,$interp->[3][1][3],$di);
    $interp->[2][3][3] = &cubefit(0.,$interp->[0][3][3],1.,$interp->[1][3][3],2.,$interp->[3][3][3],3.,$interp->[3][3][3],$di);
    $interp->[2][4][3] = &cubefit(0.,$interp->[0][4][3],1.,$interp->[1][4][3],2.,$interp->[3][4][3],3.,$interp->[3][4][3],$di);
    $interp->[2][0][4] = &cubefit(0.,$interp->[0][0][4],1.,$interp->[1][0][4],2.,$interp->[3][0][4],3.,$interp->[3][0][4],$di);
    $interp->[2][1][4] = &cubefit(0.,$interp->[0][1][4],1.,$interp->[1][1][4],2.,$interp->[3][1][4],3.,$interp->[3][1][4],$di);
    $interp->[2][3][4] = &cubefit(0.,$interp->[0][3][4],1.,$interp->[1][3][4],2.,$interp->[3][3][4],3.,$interp->[3][3][4],$di);
    $interp->[2][4][4] = &cubefit(0.,$interp->[0][4][4],1.,$interp->[1][4][4],2.,$interp->[3][4][4],3.,$interp->[3][4][4],$di);
    ##### cubically interpolate along the y direction #######
    $interp->[2][2][0] = &cubefit(0.,$interp->[2][0][0],1.,$interp->[2][1][0],2.,$interp->[2][3][0],3.,$interp->[2][4][0],$dj);
    $interp->[2][2][1] = &cubefit(0.,$interp->[2][0][1],1.,$interp->[2][1][1],2.,$interp->[2][3][1],3.,$interp->[2][4][1],$dj);
    $interp->[2][2][3] = &cubefit(0.,$interp->[2][0][3],1.,$interp->[2][1][3],2.,$interp->[2][3][3],3.,$interp->[2][4][3],$dj);
    $interp->[2][2][4] = &cubefit(0.,$interp->[2][0][4],1.,$interp->[2][1][4],2.,$interp->[2][3][4],3.,$interp->[2][4][4],$dj);
    ##### quadratically interpolate along the z direction #######
    $interp->[2][2][2] = &cubefit(0.,$interp->[2][2][0],1.,$interp->[2][2][1],2.,$interp->[2][2][3],3.,$interp->[2][2][4],$dk);

    $val = $interp->[2][2][2];
    return($val);
}