/*
DRSSTC simulator
By Antonio Carlos M. de Queiroz acmq@coe.ufrj.br
Version 1.0 - 6/6/2005
Version 1.0a - 7/6/2005 Internal steps corrected
Version 1.1 - 22/6/2005 Added vc1(0) to simulate regular Tesla coils
Version 1.1a - 26/6/2006 X0 directly calculated, G21 used (no difference)
*/

/*
This program simulates the structure:
    +v1-      k12         +v3-
 +---C1--R1--+   +--R2--+--C3--+
+|         + |   | +    |  +   |
vin       i1 L1 L2 i2   C2 v2  R3
-|         - |   | -    |  -   |
 +-----------+   +------+------+
*/

#define version "1.1a - 26/6/2005"
#include <stdio.h>
#include <conio.h>
#include <string.h>
#include <stdlib.h>
#include <ctype.h>
#include <math.h>

#define vmax 180
#define vc10 0
#define C1 5e-9
#define R1 0.0
#define L1 89.34947e-6
#define L2 28.2e-3
#define k12 0.230556
#define R2 0.0
#define C2 15e-12
#define C3 0.1e-12
#define R3 220000.0
#define fin 246420.28
#define tmax 50e-6
#define points 600
#define steps 10
#define filename "drsstc.dat"
/* Scaling factors */
#define mvc1 10
#define mvc2 1
#define mvc3 1
#define mil1 1000
#define mil2 10000

typedef double matrix [6][6];
typedef double vector [6];
int i,j,k,ii;
double G11,G12,G21,G22,M,d,dt,t;
matrix A,MT;
vector B,X0,X1,Mi={0,mvc1,mvc2,mvc3,mil1,mil2};
FILE* output;

/* Routine that inverts a matrix in place */
void Invert(int n, matrix M)
{
  int i,k,j;
  double d,e;

  for (k=1; k<=n; k++) {
    d=M[k][k]; /* May need pivotal condensation, omitted */
    if (fabs(d)<1e-12) {
      printf("Something wrong (line %d). Inversion impossible\n",k);
      exit(0);
    }
    M[k][k]=1/d;
    for (j=1; j<=n; j++)
      if (j!=k) M[k][j]=-M[k][j]/d;
    for (i=1; i<=n; i++) {
      if (i!=k) {
        e=M[i][k];
        M[i][k]=M[i][k]/d;
        for (j=1; j<=n; j++)
          if (j!=k) M[i][j]=M[i][j]+M[k][j]*e;
      }
    }
  }
}

double fvin(double t)
{
  return vmax*sin(2*M_PI*fin*t);
}

int main(void)
{
  clrscr();
  printf("DRSSTC simulator\n");
  printf("By Antonio Carlos M. de Queiroz - acmq@ufrj.br\n");
  printf("Version %s\n",version);
  printf("Element values: \n");
  printf("vmax=%g\n",vmax);
  printf("vc1(0)=%g\n",vc10);
  printf("C1=%g\n",C1);
  printf("R1=%g\n",R1);
  printf("L1=%g\n",L1);
  printf("C2=%g\n",C2);
  printf("R2=%g\n",R2);
  printf("L2=%g\n",L2);
  printf("k12=%g\n",k12);
  printf("C3=%g\n",C3);
  printf("R3=%g\n",R3);
  printf("fin=%g\n",fin);
  printf("tmax=%g\n",tmax);
  printf("points=%d\n",points);
  printf("internal steps=%d\n",steps);
  if (points*steps/(tmax*fin)<50) {
    printf("Not enough steps per cycle\n");
    exit(0);
  }
  /* Inverts [L] matrix */
  M=k12*sqrt(L1*L2);
  d=L1*L2-M*M;
  G11=L2/d;
  G22=L1/d;
  G12=-M/d;
  G21=G12;
  /* Assembles the state equations dx/dt=[A]*x+B*vin */
  A[1][1]=0; A[1][2]=0; A[1][3]=0; A[1][4]=1/C1; A[1][5]=0;
  A[2][1]=0; A[2][2]=-1/(R3*C2); A[2][3]=1/(R3*C2); A[2][4]=0; A[2][5]=-1/C2;
  A[3][1]=0; A[3][2]=1/(R3*C3); A[3][3]=-1/(R3*C3); A[3][4]=0; A[3][5]=0;
  A[4][1]=-G11; A[4][2]=G12; A[4][3]=0; A[4][4]=-G11*R1; A[4][5]=-G12*R2;
  A[5][1]=-G21; A[5][2]=G22; A[5][3]=0; A[5][4]=-G21*R1; A[5][5]=-G22*R2;
  B[1]=0; B[2]=0; B[3]=0; B[4]=G11; B[5]=G21;
  /*
  printf("State equations [A], B:\n");
  for (i=1; i<6; i++) {
    for (j=1; j<6; j++)
      printf("%12g ",A[i][j]);
    printf(",%12g\n",B[i]);
  }
  */
  dt=tmax/((points-1)*steps);
  /* Generates MT=([I]-(dt/2)*[A])^(-1) */
  for (i=1; i<6; i++)
    for (j=1; j<6; j++)
      if (i==j) MT[i][j]=1-(dt/2)*A[i][j];
      else MT[i][j]=-(dt/2)*A[i][j];
  Invert(5,MT);
  /* Sets the initial state */
  X0[1]=vc10;
  for (i=2; i<6; i++) X0[i]=0;
  t=0;
  /* Calculation */
  output=fopen(filename,"w");
  fprintf(output,"%12g",t);
  for (j=1; j<6; j++) fprintf(output," %12g",X0[j]*Mi[j]);
  fprintf(output,"\n");
  for (i=2; i<=points; i++) {
    for (ii=1; ii<=steps; ii++) { /* Internal steps */
      /* X1=X0+(dt/2)*[A]*X0+(dt/2)*B*(vin(t0)+vin(t0+dt)) */
      for (j=1; j<6; j++) {
        X1[j]=X0[j];
        for (k=1; k<6; k++) X1[j]+=(dt/2)*(A[j][k]*X0[k]);
        X1[j]+=(dt/2)*B[j]*(fvin(t)+fvin(t+dt));
      }
      /* X(t0+dt)=[MT]*X1 */
      for (j=1; j<6; j++) {
        X0[j]=0;
        for (k=1; k<6; k++) X0[j]+=MT[j][k]*X1[k];
      }
      t+=dt;
    }
    fprintf(output,"%12g",t);
    for (j=1; j<6; j++) fprintf(output," %12g",X0[j]*Mi[j]);
    fprintf(output,"\n");
  }
  fclose(output);
  printf("Output saved in file %s\n",filename);
  /* Order of the variables: v1, v2, v3, i1, i2 */
  return(1);
}