// Simulation of a three-phase fully-controlled bridge
// rectifier circuit, fed by 3-ph 4-wire input.
// This program obtainss the harmonic analysis. 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 results
// of hamonic analysis.

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

void periodicResponse(void);
void harmonicAnalysis(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;
double vSource,vRphase,vYphase,vBphase,outVolt,vCommute;
double OverLapangle,VoAvg,VoRMS,ILoadAvg,ILoadRMS;
double ILineAvg,ILineRMS,THD,RFVolt,RFCur;
int commute,NumKount,modeSw,respMode,pairNumber,fangDeg;
double anLoadCur[26],bnLoadCur[26],cnLoadCur[26];
double anLineCur[26],bnLineCur[26],cnLineCur[26];
double anOutVolt[26],bnOutVolt[26],cnOutVolt[26];

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; respMode=0;
 NumKount=0;  pairNumber=0;
 // obtain the transient response.
 periodicResponse();
 harmonicAnalysis();
 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.

 do
	{
	 startVal=cur_load;
	 OneCycle();
	} while (fabs(startVal-cur_load)>0.0005);
 // Create a file for entering transient response

}



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;
		 break;

		 case 2:
		  //SCRs 1 and 2 ON
		  vSource=vRphase-vBphase;
		  vCommute=0.0;
		 break;

		 case 3:
		  //SCRs 2 and 3 ON
		  vSource=vYphase-vBphase;
		  vCommute=0.0;
		 break;

		 case 4:
		  //SCRs 3 and 4 ON
		  vSource=vYphase-vRphase;
		  vCommute=0.0;
		 break;

		 case 5:
		  //SCRs 4 and 5 ON
		  vSource=vBphase-vRphase;
		  vCommute=0.0;
		 break;

		 case 6:
		  //SCRs 5 and 6 ON
		  vSource=vBphase-vYphase;
		  vCommute=0.0;
		 break;
		}
	}

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

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

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

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

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

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

	 case 2:
		if (commute==0)
		{
		 curRphase=cur_line;
		}
	  else
		{
		 curRphase=cur_load;
		}
	 break;

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

	 case 4:
		if (commute==0)
		{
		 curRphase=-cur_line;
		}
	  else
		{
		 curRphase=-cur_Overlap;
		}
	 break;

	 case 5:
	  if (commute==0)
		{
		 curRphase=-cur_load;
		}
	  else
		{
		 curRphase=-cur_load;
		}
	 break;

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

}


void harmonicAnalysis(void)
{
 int m;
 double theta,dblm;
 int knt;
 double dknt;
 div_t qr;
 dknt=0.0;

 for (m=0;m<26;m++)
	{
	 anLoadCur[m]=0.0; bnLoadCur[m]=0.0; cnLoadCur[m]=0.0;
	 anLineCur[m]=0.0; bnLineCur[m]=0.0; cnLineCur[m]=0.0;
	 anOutVolt[m]=0.0; bnOutVolt[m]=0.0; cnOutVolt[m]=0.0;
	}


 OverLapangle=0.0;
 ILoadAvg=0.0; ILoadRMS=0.0; RFCur=0.0;
 ILineAvg=0.0; ILineRMS=0.0; THD=0.0;
 VoAvg=0.0;  VoRMS=0.0; RFVolt=0.0;

 for (knt=0;knt<720;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();

	 if (commute==1) OverLapangle=OverLapangle+step;

	 dblm=0.0;
	 for (m=0;m<13;m++)
		{
		 theta=2.0*dblm*cycleAng;
		 while (theta>(2.0*pi)) theta=theta-2.0*pi;
		 anLoadCur[2*m]=anLoadCur[2*m]+cur_load*cos(theta)*step;
		 bnLoadCur[2*m]=bnLoadCur[2*m]+cur_load*sin(theta)*step;
		 anOutVolt[2*m]=anOutVolt[2*m]+outVolt*cos(theta)*step;
		 bnOutVolt[2*m]=bnOutVolt[2*m]+outVolt*sin(theta)*step;
		 theta=theta+cycleAng;
		 if (theta>(2.0*pi)) theta=theta-2.0*pi;
		 anLineCur[2*m+1]=anLineCur[2*m+1]+curRphase*cos(theta)*step;
		 bnLineCur[2*m+1]=bnLineCur[2*m+1]+curRphase*sin(theta)*step;
		 dblm=dblm+1.0;
		}

	 ILoadRMS=ILoadRMS+cur_load*cur_load*step;
	 VoRMS=VoRMS+outVolt*outVolt*step;
	 ILineAvg=ILineAvg+fabs(curRphase)*step;
	 ILineRMS=ILineRMS+curRphase*curRphase*step;
	 dknt+=1.0;
	}

 for (m=0;m<13;m++)
	{
	  anLoadCur[2*m]=2.0/pi*anLoadCur[2*m];
	  bnLoadCur[2*m]=2.0/pi*bnLoadCur[2*m];
	  anOutVolt[2*m]=2.0/pi*anOutVolt[2*m];
	  bnOutVolt[2*m]=2.0/pi*bnOutVolt[2*m];
	  anLineCur[2*m+1]=2.0/pi*anLineCur[2*m+1];
	  bnLineCur[2*m+1]=2.0/pi*bnLineCur[2*m+1];
	  cnLoadCur[2*m]=anLoadCur[2*m]*anLoadCur[2*m]
								 + bnLoadCur[2*m]*bnLoadCur[2*m];
	  cnLoadCur[2*m]=sqrt(cnLoadCur[2*m]);
	  cnLineCur[2*m+1]=anLineCur[2*m+1]*anLineCur[2*m+1]
								+ bnLineCur[2*m+1]*bnLineCur[2*m+1];
	  cnLineCur[2*m+1]=sqrt(cnLineCur[2*m+1]);
	  cnOutVolt[2*m]=anOutVolt[2*m]*anOutVolt[2*m]
								+ bnOutVolt[2*m]*bnOutVolt[2*m];
	  cnOutVolt[2*m]=sqrt(cnOutVolt[2*m]);
	}
 VoRMS=sqrt(VoRMS/pi);
 ILoadRMS=sqrt(ILoadRMS/pi);
 ILineAvg= ILineAvg/pi;
 ILineRMS= sqrt(ILineRMS/pi);
 ILoadAvg=anLoadCur[0]/2.0;
 VoAvg=anOutVolt[0]/2.0;
 THD=sqrt(ILineRMS*ILineRMS-cnLineCur[0]*cnLineCur[0]/2.0);
 RFVolt=sqrt(VoRMS*VoRMS-VoAvg*VoAvg);
 RFCur=sqrt(ILoadRMS*ILoadRMS-ILoadAvg*ILoadAvg);

 fnew=fopen("harm3ph.csv","w+r");
 kreateRespFile();
 produceEntries();
 fclose(fnew);
}

void kreateRespFile(void)
{
 int n1,n2;
 tstr="harmNo,LoadCur,OutVolt,harmNo,LineCur";
 p=strtok(tstr,",");
 fprintf(fnew,p); fprintf(fnew,",");
 n2=5;
 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 produceEntries(void)
{
 int m;
 for(m=0;m<13;m++)
	{
	 itoa(2*m, p, 10);
	 fprintf(fnew,p);
	 fprintf(fnew,",");

	 gcvt(cnLoadCur[2*m],5,p);
	 fprintf(fnew,p);
	 fprintf(fnew,",");

	 gcvt(cnOutVolt[2*m],5,p);
	 fprintf(fnew,p);
	 fprintf(fnew,",");

	 itoa(2*m+1, p, 10);
	 fprintf(fnew,p);
	 fprintf(fnew,",");

	 gcvt(cnLineCur[2*m+1],5,p);
	 fprintf(fnew,p);
	 fprintf(fnew,"\n");
	}

 fprintf(fnew,"\n");
 p="LoadIAvg";
 fprintf(fnew,p);
 fprintf(fnew,",");
 gcvt(ILoadAvg,5,p);
 fprintf(fnew,p);
 fprintf(fnew,"\n");
 p="LoadIRMS";
 fprintf(fnew,p);
 fprintf(fnew,",");
 gcvt(ILoadRMS,5,p);
 fprintf(fnew,p);
 fprintf(fnew,"\n");
 p="RFCur";
 fprintf(fnew,p);
 fprintf(fnew,",");
 gcvt(RFCur,5,p);
 fprintf(fnew,p);
 fprintf(fnew,"\n");
 fprintf(fnew,"\n");


 p="VoAvg";
 fprintf(fnew,p);
 fprintf(fnew,",");
 gcvt(VoAvg,5,p);
 fprintf(fnew,p);
 fprintf(fnew,"\n");
 p="VoRMS";
 fprintf(fnew,p);
 fprintf(fnew,",");
 gcvt(VoRMS,5,p);
 fprintf(fnew,p);
 fprintf(fnew,"\n");
 p="RFVolt";
 fprintf(fnew,p);
 fprintf(fnew,",");
 gcvt(RFVolt,5,p);
 fprintf(fnew,p);
 fprintf(fnew,"\n");
 fprintf(fnew,"\n");

 p="ILineAvg";
 fprintf(fnew,p);
 fprintf(fnew,",");
 gcvt(ILineAvg,5,p);
 fprintf(fnew,p);
 fprintf(fnew,"\n");
 p="ILineRMS";
 fprintf(fnew,p);
 fprintf(fnew,",");
 gcvt(ILineRMS,5,p);
 fprintf(fnew,p);
 fprintf(fnew,"\n");
 p="THD";
 fprintf(fnew,p);
 fprintf(fnew,",");
 gcvt(THD,5,p);
 fprintf(fnew,p);
 fprintf(fnew,"\n");
 fprintf(fnew,"\n");

 p="OverLapAngle";
 fprintf(fnew,p);
 fprintf(fnew,",");
 if (OverLapangle>0.0001) OverLapangle=OverLapangle/3.0+step/2.0;
 gcvt(OverLapangle,5,p);
 fprintf(fnew,p);
 fprintf(fnew,"\n");
}