Commit 968c7f21 authored by Berrill, Mark's avatar Berrill, Mark
Browse files

Syncing changes with full application

parent 58f0a5b8

ASE_medium.dat

(5.4 MiB)

File changed.

No diff preview for this file type.

View original file
View changed file

ASE_small.dat

(2.4 MiB)

File changed.

No diff preview for this file type.

View original file
View changed file

seed_medium.dat

(8.29 MiB)

File changed.

No diff preview for this file type.

View original file
View changed file

seed_small.dat

(3.71 MiB)

File changed.

No diff preview for this file type.

View original file
View changed file

src/AtomicModel/interp.cpp

+21 −17
Original line number Diff line number Diff line 
#include <algorithm>

#include <cstring>

#include <stdexcept>

#include <string.h>



#include "AtomicModel/interp.h"





static auto& perr = std::cerr;





// Subroutine to perform linear interpolation

double interp::interp_linear( size_t N, const double *xi, const double *yi, double x )

{
@@ -140,7 +143,8 @@ static inline TYPE pchip( size_t N, const TYPE *xi, const TYPE *yi, TYPE x ) 
    }

    // Perform the interpolation

    double dx2 = dx * dx;

    double f = f1 + dx2 * ( 2 * dx - 3 ) * ( f1 - f2 ) + dx * g1 - dx2 * ( g1 + ( 1 - dx ) * ( g1 + g2 ) );

    double f =

        f1 + dx2 * ( 2 * dx - 3 ) * ( f1 - f2 ) + dx * g1 - dx2 * ( g1 + ( 1 - dx ) * ( g1 + g2 ) );

    return f;

}

double interp::interp_pchip( size_t N, const double *xi, const double *yi, double x )
@@ -158,28 +162,28 @@ double interp::calc_width( size_t n, const double *x, const double *y ) 
{

    // Check the inputs

    if ( n < 2 ) {

        std::cerr << "There must be at least two points\n";

        perr << "There must be at least two points\n";

        return -1;

    }

    for ( size_t i = 0; i < n; i++ ) {

        if ( y[i] < 0.0 ) {

            std::cerr << "Negitive values in y detected\n";

            perr << "Negitive values in y detected\n";

            return -1;

        }

    }

    for ( size_t i = 1; i < n; i++ ) {

        if ( x[i] <= x[i - 1] ) {

            std::cerr << "x must be sorted and unique\n";

            perr << "x must be sorted and unique\n";

            return -1;

        }

    }

    // Calculate the normalized cumulative sum (x may not be uniformly spaced)

    double *ys = new double[n];

    auto *ys = new double[n];

    ys[0]    = 0.0;

    for ( size_t i = 1; i < n; i++ )

        ys[i] = ys[i - 1] + ( x[i] - x[i - 1] ) * 0.5 * ( y[i] + y[i - 1] ); // ys = int(y*dx)

    if ( ys[n - 1] == 0.0 ) {

        std::cerr << "y is all zeros\n";

        perr << "y is all zeros\n";

        delete[] ys;

        return -1;

    }
@@ -208,8 +212,8 @@ double interp::bisection_coeff( 
    if ( N < 2 )

        throw std::logic_error( "Error: N<2" );

    // Copy x and r

    double *x = new double[N];

    double *r = new double[N];

    auto *x = new double[N];

    auto *r = new double[N];

    memcpy( x, x_in, N * sizeof( double ) );

    memcpy( r, r_in, N * sizeof( double ) );

    // Change the sign or r such that min(x) will be at the beginning