CFUN_IcNb3Sn_HiLumi_v1.c

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

#define EPS .000000001

static const char *error = NULL;

double maximum(double a, double b) {
    return (a > b) ? a : b;
}

double minimum(double a, double b) {
    return (a < b) ? a : b;
}

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

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

EXPORT const char* return_test_yaml() {
    return "CFUN_IcNb3Sn_HiLumi_v1";
}

EXPORT int eval(const char *func,
                              int nArgs,
                              const double **inReal,
                              const double **inImag,
                              int blockSize,
                              double *outReal,
                              double *outImag) {
								  
  double Tc0   = 16.0;             //[K]
  double Bc20  = 28.11;            //[T]
  double c0    = 174000*pow(10,6); //[-]
  double degr  = 0.03;             // [-]
  double p     = 1.52;             // [-]
  double alpha = 0.96;             // [-] 
  double tmin  = 1e-2;             // Min allowed t fraction
  double Bmin  = 1e-2;             // Min allowed B
  
  double Bc2;
  double Jc;
  double ti;
  double Bi;
  double bi;
  double c;
  int i;
  
  if (strcmp("CFUN_IcNb3Sn_HiLumi_v1", func) == 0) {
    if (nArgs != 3) {
      error = "Three arguments expected";
      return 0;
    }

    for (i = 0; i < blockSize; i++) {  	 
	
	  double T    = inReal[0][i];
	  double B    = inReal[1][i];
	  double Area = inReal[2][i];
	
	  /*Input consistency check*/
	  if(T < 0)     {error = "T is negative!"; return 0; } 
	  if(B < 0)     {error = "normB is negative!"; return 0; }
	  if(Area < 0)  {error = "Area is negative!"; return 0; }	
	  
       /* Mandatory constraint: 0.01 <= t <= 1 */
       ti=T/Tc0;
	   
	   if(ti>=1) {
		   Jc = 0;
	   }	   
	   else {
			   
		   ti=maximum(ti,tmin);
		   
		   Bc2 = Bc20 * (1-pow(ti,p)); 

		   /* Mandatory constraint: 0.01 <= b <= 1 */
		   bi=B/Bc2;
		   bi=minimum(maximum(bi,Bmin),1);       

		   /* Bordini's original fit: */
		   c = (1-degr) * c0 * pow((1-pow(ti,p)),alpha) * pow((1-pow(ti,2)),alpha);
		  Jc = c / (Bc2*bi) * sqrt(bi)*pow((1-bi),2);	
	   }

       /* Note that this equation may also be mathematically simplified to:
       * c = (1-degr) * c0/Bc20 * pow((1-pow(ti,p)),alpha-1) * pow((1-pow(ti,2)),alpha);
       * Jc = c * pow(bi,-0.5) * pow((1-bi),2)
       */

       outReal[i] = Jc*Area;
	   
	   	  /*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;

  }

  else {
    error = "Unknown function";
    return 0;
  }
}