// Simulation of a three-phase fully-controlled bridge
// rectifier circuit, fed by 3-ph 4-wire input.
// This program obtainss the periodic response. The program  needs
// only three parameters to be keyed in: the ratio of load
// inductive reactance to the load resistance, wL1/R , the
// ratio of source reactance to the load resistance to load
// resistance wL2/R and the firing angle in degrees.
// The program produces an output file containing the response
// obtained over a cycle.  The parameters stored are :
// angle ranging from 0 to 360 deg, phase R input voltage,
// phase R  current, the voltage across SCR S1 connected
// between R phase and positive rail, the voltage output from
// the bridge, the load current, the voltage across the line  or
// ac choke and the voltage across the dc link inductance.

#include <math.h>
#include<stdlib.h>
#include<stdio.h>
#include<string.h>

void periodicResponse(void);
void OneCycle(void);
void computeOneStep(void);
void kreateRespFile(void);
void produceEntries(void);

const double pi=3.1415926;
const double deg_rad=pi/180.0;
const double step=pi/720.0;
const double MAG=pi/3.0/sqrt(3.0);

double L1,L2,fangRad,cycleAng;
double cur_load,cur_line,curRphase,cur_Overlap,acIndV,dcIndV;
double vSource,vRphase,vYphase,vBphase,outVolt,vCommute,vSCR1;
int commute,NumKount,modeSw,respMode,pairNumber,fangDeg;
FILE *fnew;
char *tstr,*p;


void main(void)
{
 float  ka1,kb2,kc3;
 printf(" \n");
 printf("     Load Time Constant in radians = ");
 scanf("%f",&ka1);
 L1=(double)ka1;
 if (L1<0.05) L1=0.1;
 printf(" \n");
 printf("     Line Reactance Time Constant in radians = ");
 scanf("%f",&kb2);
 L2=(double)kb2;
 printf(" \n");
 printf("     Firing angle in degrees = ");
 scanf("%f",&kc3);
 fangRad=deg_rad*(double)kc3;
 fangDeg=(int)kc3;

 printf("  \n");


 // Initialise though not necessary for global variables

 cycleAng=0.0;
 cur_load=0.0;  cur_line=0.0;  outVolt=0.0;
 commute=0; modeSw=0;
 NumKount=0;  pairNumber=0;
 // obtain the transient response.
 periodicResponse();

 printf("Close window by clicking on x button at top right corner \n");
}

void periodicResponse(void)
{
 double startVal;

 // Calculate one cycle of response and repeat till the
 // response becomes periodic.
 respMode=0;
 do
	{
	 startVal=cur_load;
	 OneCycle();
	} while (fabs(startVal-cur_load)>0.0005);
 // Create a file for entering transient response
 fnew=fopen("ss_resp1.csv","w+r");
 kreateRespFile();
 respMode=1;
 OneCycle();
 fclose(fnew);
}

void kreateRespFile(void)
{
 int n1,n2;
 tstr="Angle,ph_RVolt,ph_RCur,vSCR1,VoutBr,LoadCur,acIndV,dcIndV";
 p=strtok(tstr,",");
 fprintf(fnew,p); fprintf(fnew,",");
 n2=8;
 for (n1=0;n1<(n2-1);n1++)
	{
	 p=strtok(NULL,",");
	 if (n1!=(n2-2))
		{
		 fprintf(fnew,p); fprintf(fnew,",");
		}
	 else
		{
		 fprintf(fnew,p); fprintf(fnew,"\n");
		}
	}
}


void OneCycle(void)
{
 int knt;
 double dknt;
 dknt=0.0;
 div_t qr;

 for (knt=0;knt<1440;knt++)
	{
	 qr=div(knt,4);
	 if ((qr.quot>=(30+fangDeg)) &&(qr.quot<(90+fangDeg)))
												 pairNumber=1; //S6,S1 ON
	 if ((qr.quot>=(90+fangDeg)) &&(qr.quot<(150+fangDeg)))
												 pairNumber=2; //S1,S2 ON
	 if ((qr.quot>=(150+fangDeg)) &&(qr.quot<(210+fangDeg)))
												 pairNumber=3; //S2,S3 ON
	 if ((qr.quot>=(210+fangDeg)) &&(qr.quot<(270+fangDeg)))
												 pairNumber=4; //S3,S4 ON
	 if (fangDeg<30)
		{
		 if ((qr.quot>=(270+fangDeg)) &&(qr.quot<(330+fangDeg)))
												 pairNumber=5; //S4,S5 ON
		 if ((qr.quot>=(330+fangDeg)) || ((qr.quot+360)<(390+fangDeg)))
												 pairNumber=6; //S5,S6 ON
		}
	 if ((fangDeg>=30) && (fangDeg<90))
		{
		 if ((qr.quot>=(270+fangDeg)) && (qr.quot<(330+fangDeg)))
												 pairNumber=5; //S5,S6 ON
       if ((qr.quot+360)<(330+fangDeg) && (pairNumber==5)) pairNumber=5;
		 if (((qr.quot+360)>=(330+fangDeg)) && ((qr.quot+360)<(390+fangDeg)))
												 pairNumber=6; //S6,S1 ON
		}
	 cycleAng=dknt*step;

	 // routine to start commutation overlap
	 // When the current through the incoming SCR equals
	 // the load current, commutation overlap ends.
	 if (modeSw!=pairNumber)
		{
		 commute=1;
		 if (L2<0.0005) commute=0;
		 modeSw=pairNumber;
		}
	 computeOneStep();
	 dknt+=1.0;
	}
}

void computeOneStep(void)
{
 double dLoadI,dLineI;

 vRphase = MAG*sin(cycleAng);
 vYphase = MAG*sin(cycleAng-2.0*pi/3.0);
 vBphase = MAG*sin(cycleAng+2.0*pi/3.0);

 if (commute==0)
	{
	 switch (pairNumber)
		{
		 case 1:
		  //SCRs 1 and 6 ON
		  vSource=vRphase-vYphase;
		  vCommute=0.0;
		  if (cur_load>0.0) vSCR1=0.0;
		  else vSCR1=vRphase;
		 break;

		 case 2:
		  //SCRs 1 and 2 ON
		  vSource=vRphase-vBphase;
		  vCommute=0.0;
		  if (cur_load>0.0) vSCR1=0.0;
		  else vSCR1=vRphase;
		 break;

		 case 3:
		  //SCRs 2 and 3 ON
		  vSource=vYphase-vBphase;
		  vCommute=0.0;
		  if (cur_load>0.0) vSCR1=vRphase-vYphase;
		  else vSCR1=vRphase;
		 break;

		 case 4:
		  //SCRs 3 and 4 ON
		  vSource=vYphase-vRphase;
		  vCommute=0.0;
		  if (cur_load>0.0) vSCR1=vRphase-vYphase;
		  else vSCR1=vRphase;
		 break;

		 case 5:
		  //SCRs 4 and 5 ON
		  vSource=vBphase-vRphase;
		  vCommute=0.0;
		  if (cur_load>0.0) vSCR1=vRphase-vBphase;
		  else vSCR1=vRphase;
		 break;

		 case 6:
		  //SCRs 5 and 6 ON
		  vSource=vBphase-vYphase;
		  vCommute=0.0;
		  if (cur_load>0.0) vSCR1=vRphase-vBphase;
		  else vSCR1=vRphase;
		 break;
		}
	}

 else
	{
	 switch (pairNumber)
		{
		 case 1:
		  // SCR 6 ON. SCR1 takes over from SCR5.
		  vSource=(vRphase+vBphase)/2.0-vYphase;
		  vCommute=vRphase - vBphase;
		  vSCR1=0.0;
		 break;

		 case 2:
		  // SCR 1 ON. SCR2 takes over from SCR6.
		  vSource=vRphase-(vBphase+vYphase)/2.0;
		  vCommute= vYphase - vBphase;
		  vSCR1=0.0;
		 break;

		 case 3:
		  // SCR 2 ON. SCR3 takes over from SCR1.
		  vSource=(vYphase+vRphase)/2.0-vBphase;
		  vCommute=vYphase - vRphase;
		  vSCR1=0.0;
		 break;

		 case 4:
		  // SCR 3 ON. SCR4 takes over from SCR2.
		  vSource=vYphase-(vBphase+vRphase)/2.0;
		  vCommute=vBphase - vRphase;
		  vSCR1=vRphase-vYphase;
		 break;

		 case 5:
		  // SCR 4 ON. SCR5 takes over from SCR3.
		  vSource=(vBphase+vYphase)/2.0-vRphase;
		  vCommute= vBphase-vYphase;
		  vSCR1=vRphase-(vYphase+vBphase)/2.0;
		 break;

		 case 6:
		  // SCR 5 ON. SCR6 takes over from SCR4.
		  vSource=vBphase-(vRphase+vYphase)/2.0;
		  vCommute=vRphase - vYphase;
		  vSCR1=vRphase-vBphase;
		 break;
		}
	}
 switch (commute)
	{
	 case 0:
	  // cur_line and cur_load are equal.
	  cur_Overlap=0.0;
	  dLoadI=(vSource-cur_load)/(L1+2.0*L2)*step;
	  cur_load=cur_load+dLoadI;
	  if (cur_load<0.00001) cur_load=0.0;
	  cur_line=cur_load;
	  if (cur_load>0.0) outVolt=cur_load+(vSource-cur_load)*L1/(L1+2.0*L2);
	  else outVolt=0.0;
	 break;

	 case 1:
	  // Current through the SCR turned on is computed.
	  // When it equals, load current, commutation is over.
	  dLoadI=(vSource-cur_load)/(L1+1.5*L2)*step;
	  cur_load=cur_load+dLoadI;
	  if (cur_load<0.00001) cur_load=0.0;
	  dLineI=vCommute/(2.0*L2)*step;
	  cur_Overlap=cur_Overlap+dLineI;
	  if (cur_Overlap>cur_load) commute=0;
	  outVolt=cur_load+(vSource-cur_load)*L1/(L1+1.5*L2);
	 break;
	}

 // calculate current in R phase.  This information can be
 // obtained from the load current and the pair in conduction
 // when there is no overlap.  If there is overlap, we need to
 // know the load current, the pair in conduction and
 // the current through the SCR being turned on.

 switch (pairNumber)
	{
	 case 1:
	  if (commute==0)
		{
		 curRphase=cur_line;
		 acIndV=(vSource-outVolt)/2.0;
		}
	  else
		{
		 curRphase=cur_Overlap;
		 acIndV=vCommute/2.0;
		}
	 break;

	 case 2:
		if (commute==0)
		{
		 curRphase=cur_line;
		 acIndV=(vSource-outVolt)/2.0;
		}
	  else
		{
		 curRphase=cur_load;
		 acIndV=(vSource-outVolt)*2.0/3.0;
		}
	 break;

	 case 3:
	  if (commute==0)
		{
		 curRphase=0.0;
		 acIndV=0.0;
		}
	  else
		{
		 curRphase=cur_load-cur_Overlap;
		 acIndV=-vCommute/2.0;
		}
	 break;

	 case 4:
		if (commute==0)
		{
		 curRphase=-cur_line;
		 acIndV=-(vSource-outVolt)/2.0;
		}
	  else
		{
		 curRphase=-cur_Overlap;
		 acIndV=-vCommute/2.0;
		}
	 break;

	 case 5:
	  if (commute==0)
		{
		 curRphase=-cur_load;
		 acIndV=-(vSource-outVolt)/2.0;
		}
	  else
		{
		 curRphase=-cur_load;
		 acIndV=-(vSource-outVolt)*2.0/3.0;
		}
	 break;

	 case 6:
	  if (commute==0)
		{
		 curRphase=0.0;
		 acIndV=0.0;
		}
	  else
		{
		 curRphase=-(cur_load-cur_Overlap);
		 acIndV=vCommute/2.0;
		}
	 break;
	}

 dcIndV=outVolt-cur_load;
 if (respMode==1) produceEntries();
}


void produceEntries(void)
{

 gcvt(cycleAng/pi*180.0,5,p);
 fprintf(fnew,p);
 fprintf(fnew,",");

 gcvt(vRphase,5,p);
 fprintf(fnew,p);
 fprintf(fnew,",");

 gcvt(curRphase,5,p);
 fprintf(fnew,p);
 fprintf(fnew,",");

 gcvt(vSCR1,5,p);
 fprintf(fnew,p);
 fprintf(fnew,",");


 gcvt(outVolt,5,p);
 fprintf(fnew,p);
 fprintf(fnew,",");

 gcvt(cur_load,5,p);
 fprintf(fnew,p);
 fprintf(fnew,",");

 gcvt(acIndV,5,p);
 fprintf(fnew,p);
 fprintf(fnew,",");

 gcvt(dcIndV,5,p);
 fprintf(fnew,p);
 fprintf(fnew,"\n");

}