CFUN_HTS_JcFit_SUPERPOWER_v1_deriv.c

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

#define MAX(x, y) (((x) > (y)) ? (x) : (y))
#define MIN(x, y) (((x) < (y)) ? (x) : (y))
#define PI 3.14159265358979323846
#define deg2rad PI/180
#define rad2deg 180/PI

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) {
 
 	
 if (nArgs != 3) {
	// Top_IN, Bnorm_IN, thetaFieldTape_IN
	error = "three arguments expected";
    return 0;
    }
	
// Fit Parameters
double n0;  
double n1;
double n2;
double a; 
double Tc0;
double Bi0_ab;
double Bi0_c;
double alpha_ab;
double p_ab;
double q_ab;
double alpha_c;
double p_c;
double q_c;
double g0;
double g1;
double g2;
double g3;
double gamma_ab;
double gamma_c;
double nu;
double thetaOff;

if(strcmp("CFUN_HTS_JcFit_SUPERPOWER_v1_deriv", func) == 0) {
	g0       = 0.039445688; //[-]
	g1       = 0.276598888;  //[-]
	g2       = 0.047699303; //[-]
	g3       = 0.055878401; //[-]
	Tc0      = 93; //[K]
	p_c      = 0.400002709; //[-]
	q_c      = 2.940671164; //[-]
	Bi0_c    = 140; //[T]
	gamma_c  = 2.999950877; //[-]
	alpha_c  = 2.99996E+12; //[A*T/m^2]
	nu       = 1.720055527; //[-]		
	n0       = 1.169541488; //[-]
	n1       = 1.599990183; //[-]
	n2       = 1.131992777; //[-]
	p_ab     = 0.825413689; //[-]
	q_ab     = 5.495925721; //[-]
	Bi0_ab   = 250; //[T]		
	a        = 2.05E-05; //[-]
	gamma_ab = 4.998961387; //[-]
	alpha_ab = 7.10415E+13; //[A*T/m^2]
	thetaOff = -2.179824999; //[deg]	
}

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

// Correction factor
double B_MIN = 0.01; //[T]

// Global index
int i; 

// Parfor loop
for (i = 0; i < blockSize; i++) {
	// Raw Input
	double Top_IN   = inReal[0][i];	//[0,100]
	double Bnorm_IN = inReal[1][i]; // 4.1
	double thetaFieldTape_IN = inReal[2][i]; // 0

	//Pre-allocation
	double Top;
	double Bnorm;
	double theta;
	double thetaMod;
	
	//Pre-processing
	Bnorm     = MAX(B_MIN,fabs(Bnorm_IN));
	Top       = Top_IN;	
    theta     = fabs(thetaFieldTape_IN - thetaOff)*deg2rad;	

	if (fmod(theta,PI)<=PI/2){ 
		thetaMod = fmod(theta,PI); 
		}
	else{ 
		thetaMod = PI-fmod(theta,PI);
		}

	//Pre-allocation
	double t_T;
	double t_n0;	
	double t_n1;	
	double Bi_ab;
	double Bi_c;
	double b_ab;
	double b_c;
	double Jc_ab;
	double Jc_c;
	bool TCrit_Flag;	
	bool babCrit_Flag;
	bool bcCrit_Flag;
	double g;
	double JcFit;
	double dt_T, dTop, dt_n0_dT, dg_dT, dt_n1_dT, dJc_ab_dt_T, derivative_Jc_ab, derivative_Jc_c, JcFit_derivative, dBi_ab_dT, dBi_c_dT, db_ab_dT, db_c_dT;
	
	// Temperature scaling factors
    t_T  = Top/Tc0;
	dt_T = 1/Tc0;

	t_n0 = 1-pow(t_T,n0);
	dt_n0_dT = - n0 * pow(t_T,n0-1) * dt_T;
	
	t_n1 = 1-pow(t_T,n1);	
	dt_n1_dT = - n1 * pow(t_T,n1-1) * dt_T;

	// Irreversibility field  
	Bi_ab = Bi0_ab*(pow(t_n1,n2)+a*t_n0);
	dBi_ab_dT = Bi0_ab*(n2* pow(t_n1,n2 -1) * dt_n1_dT +a*dt_n0_dT);

	Bi_c  = Bi0_c*t_n0;
	dBi_c_dT  = Bi0_c*dt_n0_dT;
    
	// magnetic flux density	
	b_ab  = Bnorm / Bi_ab;
	db_ab_dT = - Bnorm / pow(Bi_ab, 2) *dBi_ab_dT;
	b_c   = Bnorm / Bi_c;
	db_c_dT = - Bnorm / pow(Bi_c, 2) *dBi_c_dT;
	
	// flags	
	TCrit_Flag   = t_T  >= 1 || t_T  < 0;
	babCrit_Flag = b_ab >= 1 || b_ab < 0;		
	bcCrit_Flag  = b_c  >= 1 || b_c  < 0;
	
	if (TCrit_Flag==false){
			
		if (babCrit_Flag==false){			
			derivative_Jc_ab =	(alpha_ab/Bnorm*pow(b_ab,p_ab)) * pow((1-b_ab),q_ab) * ((n2 *pow(t_n1,n2-1)* dt_n1_dT + a* dt_n0_dT)*gamma_ab *pow((pow(t_n1,n2) + a*t_n0),gamma_ab - 1)) + (alpha_ab/Bnorm* (p_ab*pow(b_ab,p_ab-1)) * pow((1-b_ab),q_ab)+(alpha_ab/Bnorm*pow(b_ab,p_ab)) *(-q_ab)* pow((1-b_ab),q_ab-1))*db_ab_dT* pow((pow(t_n1,n2) + a*t_n0),gamma_ab);

			Jc_ab =	(alpha_ab/Bnorm*pow(b_ab,p_ab)) * pow((1-b_ab),q_ab) * pow((pow(t_n1,n2) + a*t_n0),gamma_ab);
			}
		else{
			// field over Bi
			Jc_ab = 0;
			derivative_Jc_ab = 0;
			}
		
		if (bcCrit_Flag==false){
			// critical current density c-plane
			derivative_Jc_c = (1/Bnorm)*alpha_c*pow(b_c,p_c)*pow((1-b_c),q_c)* gamma_c * pow((1-pow(t_T,n0)),gamma_c-1) * dt_T * n0 *(-pow(t_T,n0-1)) + (1/Bnorm)*alpha_c*(p_c*pow(b_c,p_c-1)*pow((1-b_c),q_c) - q_c*pow(b_c,p_c)*pow((1-b_c),q_c-1))*pow((1-pow(t_T,n0)),gamma_c)*db_c_dT;

			Jc_c = (1/Bnorm)*alpha_c*pow(b_c,p_c)*pow((1-b_c),q_c)*pow((1-pow(t_T,n0)),gamma_c);
			}
		else{
			// field over Bi
			Jc_c = 0;
			derivative_Jc_c = 0;
			}	
        
		// anisotropy factor
		g = g0 + g1*exp(-g2*Bnorm*exp(g3*Top));
		dg_dT = -g1 * g2 * g3 * Bnorm *exp(g3 * Top)* exp(-g2 * Bnorm * exp(g3 * Top));
		//g1*exp(-g2*Bnorm*exp(g3*Top));

		if(Jc_ab > Jc_c){
			//JcFit_derivative = derivative_Jc_c + (derivative_Jc_ab - derivative_Jc_c) /(1+pow(((PI/2-thetaMod)/g),nu)) + (Jc_ab-Jc_c)*dg_dT*pow((PI/2-thetaMod), nu)*(-nu)*pow(1/g,nu+1)*(-1)/pow((1+pow((PI/2-thetaMod)/g, nu)),2);

			JcFit_derivative = derivative_Jc_c + (derivative_Jc_ab - derivative_Jc_c) /(1+pow(((PI/2-thetaMod)/g),nu)) + (Jc_ab-Jc_c)*dg_dT/pow((1+pow((PI/2-theta)/g, nu)), 2) * (-1)* nu * pow((PI/2-theta)/g, nu-1) * (-(PI/2-theta))/pow(g,2);
		}
		else if (Jc_ab < Jc_c){
			JcFit_derivative = derivative_Jc_ab;
		}
		else{
		// temperature over Tc0
		JcFit_derivative=0;
		}
	}
		else{
		// temperature over Tc0
		JcFit_derivative=0;
		}

	//outReal[i] = JcFit;
	outReal[i] =JcFit_derivative;
	
	//Consistency check: output
	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;
}