#include <math.h>
#include <float.h>
#include <stdlib.h>
#include <string.h>
#ifdef _MSC_VER
#define EXPORT __declspec(dllexport)
#else
#define EXPORT
#endif

static const char *error = NULL;

EXPORT int init(const char *str) {
  return 1;
}

EXPORT const char * getLastError() {
  return error;
}


EXPORT int eval(const char *func,
                              int nArgs,
                              const double **inReal,
                              const double **inImag,
                              int blockSize,
                              double *outReal,
                              double *outImag) {
    
    int i;
    double a_0 = 21.945;

    double b_0 = -131.22;
    double b_1 = 37.64 * 2;
    double b_2 = -1.13 * 3;

    double c_0 = -723.04;
    double c_1 = 13.48 * 2;
    double c_2 = -0.084 * 3;

    double d_0 = 0.069;

    double threshold_0 = 0;
    double threshold_1 = 4;
    double threshold_2 = 19.6;
    double threshold_3 = 52;
    double threshold_4 = 214;

    if (strcmp("CFUN_hN_v1_deriv", func) == 0) {
        if (nArgs != 1) {
        error = "One argument expected";
        return 0;}
        }
    //It returns heat transfer coefficient of liquid Nitrogen bath in units [W/(K*m^2)]
    // It is fitting II of the following paper
    // Source: Jin, T., Hong, Jp., Zheng, H. et al. 
    // Measurement of boiling heat transfer coefficient in liquid nitrogen bath by inverse heat conduction method. 
    // J. Zhejiang Univ. Sci. A 10, 691–696 (2009). 
    // https://doi.org/10.1631/jzus.A0820540
    for (i = 0; i < blockSize; i++) {
      double T_diff        = inReal[0][i]; // Difference between the temperature of the object and the liquid nitrogen bath

      if(T_diff <= threshold_0) 
        outReal[i] = 0.0; 
      else if (T_diff < threshold_1)
        outReal[i] = a_0;
      else if (T_diff < threshold_2)
        outReal[i] = b_0 + b_1 * T_diff + b_2 * pow(T_diff, 2);
      else if (T_diff < threshold_3)
        outReal[i] = c_0 + c_1 * T_diff + c_2 * pow(T_diff, 2);
      else if (T_diff < threshold_4)
        outReal[i] = d_0;
      else //if above threshold_4, the value stays constant, this happens when the object is above 77K + 214K = 291K
        outReal[i] = 0;
      
        /*Output consistency check*/
      if(outReal[i]!=outReal[i]){
        error = "Output is nan"; 
        return 0;	  
        }
      if (fabs(outReal[i])>DBL_MAX){
        error = "Output is inf"; 
        return 0;	  
        }
    }
    return 1;
}