CFUN_QHCircuit_v3.c

#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 const char* return_test_yaml() {
    return "CFUN_QHCircuit_v3";
}

EXPORT int eval(const char *func,
                              int nArgs,
                              const double **inReal,
                              const double **inImag,
                              int blockSize,
                              double *outReal,
                              double *outImag) {
  int i;
  
  if (strcmp("CFUN_QHCircuit_v3", func) == 0) {
    if (nArgs != 11) {
      error = "Eleven arguments expected";
      return 0;
    }

    for (i = 0; i < blockSize; i++) {
	  /*Input parameters read*/
	  double t    = inReal[0][i];    // Time [s]
	  double rho_SS = inReal[1][i];  // Resistivity of stainless steel [Ohm m]
	  double t_on = inReal[2][i];    // Time when the QH system is triggered (t=0 of quench detection) [s]
	  double U_0    = inReal[3][i];  // Charging voltage of the QH capacitor bank [V]
	  double C      = inReal[4][i];  // Capacitance of the QH capacitor bank [F]
	  double R_warm = inReal[5][i];  // Additional resistance in the room temperature part of the QH circuit [Ohm]
	  double w_SS   = inReal[6][i];  // Width of the QH strip [m]
	  double h_SS   = inReal[7][i];  // Thickness of the stainless steel portion of the QH strip [m]
	  double l_SS   = inReal[8][i];  // Length of the total stainless steel portion of the QH strip [m]
	  double L   = inReal[9][i];  	 // Inductance of the QH circuit [H] 
	  int mode      = (int) inReal[10][i];

	  
	  /*Local variables needed for computation*/
	  double i_0, i_SS, tau, S_SS, R_SS,R_total,alpha,omega_0,omega_d,lambda1,lambda2,s;
	  
	  /*Input consistency check*/
	  if(t<0)     {error = "Time is negative!";        return 0; }
	  if(t_on<0)  {error = "t_on should be positive!"; return 0; }
	  if(C<0)     {error = "Capacitance is negative!"; return 0; }
	  if(R_warm<0){error = "R_warm is negative!";      return 0; }
	  if(w_SS<0)  {error = "w_SS is negative!";        return 0; }
	  if(h_SS<0)  {error = "h_SS is negative!";        return 0; }
	  if(l_SS<0)  {error = "l_SS is negative!";        return 0; }
	  if(L<0)  	  {error = "L is negative!";           return 0; }
	  
	  if(t >= t_on) {
		  t    = t - t_on;
		  S_SS = (w_SS * h_SS);

		  R_SS    = rho_SS * l_SS / S_SS;

         R_total = R_warm + R_SS;
         alpha = R_total / (2.0 * L);
         omega_0 = 1.0 / sqrt(L * C);


        if (alpha*alpha < omega_0*omega_0) {
            // Underdamped
            omega_d = sqrt(omega_0 * omega_0 - alpha * alpha);
            i_SS = (U_0 / (omega_d * L)) * exp(-alpha * t) * sin(omega_d * t);
        } else if (alpha == omega_0) {
            // Critically damped
            i_SS = (U_0 / L) * t * exp(-alpha * t);
        } else {
            // Overdamped
            s = sqrt(alpha * alpha - omega_0 * omega_0);
             lambda1 = alpha + s;
             lambda2 = alpha - s;
            i_SS = (U_0 / (L * (lambda2 - lambda1))) * (exp(-lambda1 * t) - exp(-lambda2 * t));
        }



		  switch(mode) {
			case 1: // Power
				outReal[i] = pow(i_SS / S_SS, 2) * rho_SS;
				break;
			case 2: // Current
				outReal[i] = i_SS;
				break;
			case 3: // Resistance
				outReal[i] = R_SS;
				break;
			default:
				outReal[i] = 0.0;
				break;
		  }
	  }
	  else {
		  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;
  }

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