-
Leon Teichrob authoredLeon Teichrob authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
CFUN_HTS_JcFit_SUPERPOWER_v1_deriv.c 5.67 KiB
#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;
}