/*
TRSSTC simulator
By Antonio Carlos M. de Queiroz acmq@coe.ufrj.br
Version 1.0 - 22/6/2005
Version 1.0a - 26/6/2005 X0 calculated directly, G21 used (no difference)
*/

/*
This program simulates the structure:
    +v1-      k12         -i3+
 +---C1--R1--+   +--R2--+--L3--R3--+--R4--+
+|         + |   | +    |  +       |  +   |  +
vin       i1 L1 L2 i2   C2 v2      C3 v3  C4 v4
-|         - |   | -    |  -       |  -   |  -
 +-----------+   +------+----------+------+
*/
/* To reduce to a DRSSTC:
Split L2 into L2 and L3
Multiply k12 by sqrt((L2+L3)/L2)
Split R2 into R2 and R3
Make C2 small
*/

#define version "1.0 - 22/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 C2 0.01e-12
#define R2 0.0
#define L2 14.1e-3
#define k12 0.326055
#define C3 15e-12
#define R3 0.0
#define L3 14.1e-3
#define C4 0.1e-12
#define R4 220000.0
#define fin 246420.28
#define tmax 50e-6
#define points 600
#define steps 10
#define filename "trsstc.dat"
/* Scaling factors */
#define mvc1 10
#define mvc2 1
#define mvc3 1
#define mvc4 1
#define mil1 1000
#define mil2 1000
#define mil3 10000

typedef double matrix [8][8];
typedef double vector [8];
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,mvc4,mil1,mil2,mil3};
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("TRSSTC 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("L3=%g\n",L3);
  printf("C4=%g\n",C4);
  printf("R4=%g\n",R4);
  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]=0; A[1][5]=1/C1; A[1][6]=0; A[1][7]=0;
  A[2][1]=0; A[2][2]=0; A[2][3]=0; A[2][4]=0; A[2][5]=0; A[2][6]=-1/C2; A[2][7]=1/C2;
  A[3][1]=0; A[3][2]=0; A[3][3]=-1/(R4*C3); A[3][4]=1/(R4*C3); A[3][5]=0; A[3][6]=0; A[3][7]=-1/C3;
  A[4][1]=0; A[4][2]=0; A[4][3]=1/(R4*C4); A[4][4]=-1/(R4*C4); A[4][5]=0; A[4][6]=0; A[4][7]=0;
  A[5][1]=-G11; A[5][2]=G12; A[5][3]=0; A[5][4]=0; A[5][5]=-G11*R1; A[5][6]=-G12*R2; A[5][7]=0;
  A[6][1]=-G21; A[6][2]=G22; A[6][3]=0; A[6][4]=0; A[6][5]=-G21*R1; A[6][6]=-G22*R2; A[6][7]=0;
  A[7][1]=0; A[7][2]=-1/L3; A[7][3]=1/L3; A[7][4]=0; A[7][5]=0; A[7][6]=0; A[7][7]=-R3/L3;
  B[1]=0; B[2]=0; B[3]=0; B[4]=0; B[5]=G11; B[6]=G21; B[7]=0;
  /*
  printf("State equations [A], B:\n");
  for (i=1; i<8; i++) {
    for (j=1; j<8; 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<8; i++)
    for (j=1; j<8; j++)
      if (i==j) MT[i][j]=1-(dt/2)*A[i][j];
      else MT[i][j]=-(dt/2)*A[i][j];
  Invert(7,MT);
  /* Sets the initial state */
  X0[1]=vc10;
  for (i=2; i<8; i++) X0[i]=0;
  t=0;
  /* Calculation */
  output=fopen(filename,"w");
  fprintf(output,"%12g",t);
  for (j=1; j<8; 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<8; j++) {
        X1[j]=X0[j];
        for (k=1; k<8; 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<8; j++) {
        X0[j]=0;
        for (k=1; k<8; k++) X0[j]+=MT[j][k]*X1[k];
      }
      t+=dt;
    }
    fprintf(output,"%12g",t);
    for (j=1; j<8; 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, v4, i1, i2, i3 */
  return(1);
}