// 潮流.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include <iostream>
#include <iomanip>
using namespace std;

void main()

{	int i,j,k=0;
	double x[12]={0.1,0,0.12,0.1,0,0.08,0,0,0,0.12,0.08,0},
		   y[12]={0.4,0.27272727,0.5,0.4,0,0.4,0.27272727,0,0,0.5,0.4,0},
		   z[12]={0.01528,-0.11/0.3,0.0192,0.01528,0,0.01413,1.0/3.0,0,0,0.0192,0.01413,0};
	double m[4][4];//節(jié)點導(dǎo)納矩陣實部
	double n[4][4];//節(jié)點導(dǎo)納矩陣虛部
	for(i=0;i<4;i++)
		for(j=0;j<4;j++)
		{
			m[i][j]=0;
			n[i][j]=0;
		}
	for(i=0;i<=3;i++)
	{
		for(j=0;j<=3;j++)
		{
		    if(i>j)
				if(x[i*3+j]*x[i*3+j]+y[i*3+j]*y[i*3+j]==0)
				    {m[i][j]=0;n[i][j]=0;}
				else
					{m[i][j]=-x[i*3+j]/(x[i*3+j]*x[i*3+j]+y[i*3+j]*y[i*3+j]);
				     n[i][j]=y[i*3+j]/(x[i*3+j]*x[i*3+j]+y[i*3+j]*y[i*3+j]);}
			else if(i<j)
				if(x[i*3+j-1]*x[i*3+j-1]+y[i*3+j-1]*y[i*3+j-1]==0)
					{m[i][j]=0;n[i][j]=0;}
			    else
				{m[i][j]=-x[i*3+j-1]/(x[i*3+j-1]*x[i*3+j-1]+y[i*3+j-1]*y[i*3+j-1]);
				 n[i][j]=y[i*3+j-1]/(x[i*3+j-1]*x[i*3+j-1]+y[i*3+j-1]*y[i*3+j-1]);}
		    else
			{
				for(k=0;k<3;k++)
				{
					if((x[i*3+k]*x[i*3+k]+y[i*3+k]*y[i*3+k])==0)
					{m[i][j]+=0;n[i][j]+=0;}
					else
					{m[i][j]+=x[i*3+k]/(x[i*3+k]*x[i*3+k]+y[i*3+k]*y[i*3+k]);
					n[i][j]+=-y[i*3+k]/(x[i*3+k]*x[i*3+k]+y[i*3+k]*y[i*3+k])+z[i*3+k];}
				}}}}
	cout<<"節(jié)點導(dǎo)納矩陣"<<endl;
	for(i=0;i<=3;i++)
	{for(j=0;j<=3;j++)
	 cout<<setw(12)<<setprecision(6)<<m[i][j]<<"+"<<setprecision(7)<<n[i][j];}
double ps[4]={-0.3,-0.55,0.5},
       qs[4]={-0.18,-0.13},
	   e[4]={1,1,1.1,1.05},
	   f[4]={0,0,0,0},
	   v3s=1.1,v4s=1.05,
	   p[3],q[2],vv;
double X=0,Y=0,C=0,D=0;

	for(i=0;i<=2;i++)
 {  X=0;
    Y=0;
	for(j=0;j<=3;j++)
	{X+=m[i][j]*e[j]-n[i][j]*f[j];  
     Y+=m[i][j]*f[j]+n[i][j]*e[j];
	}
	p[i]=ps[i]-(e[i]*X+f[i]*Y);
  }
for(i=0;i<=1;i++)
  { C=0;
    D=0;
	for(j=0;j<=3;j++)
	{C+=m[i][j]*e[j]-n[i][j]*f[j];  
     D+=m[i][j]*f[j]+n[i][j]*e[j];}
	q[i]=qs[i]-(f[i]*C-e[i]*D);
  }
vv=v3s*v3s-(e[2]*e[2]+f[2]*f[2]);
cout<<"ΔP ΔQ ΔV的值"<<endl;
for(i=0;i<=2;i++)
cout<<p[i]<<" ";
cout<<endl;
for(i=0;i<=1;i++)
cout<<q[i]<<" ";
cout<<endl;
cout<<vv<<endl;
double J[6][6];
for(i=0;i<6;i++)      //求雅可比行列式
	{for(int j=0;j<6;j++)
	J[i][j]=0;}
for(i=0;i<6;i++)
	{if(i%2==0)
	{for(int j=0;j<6;j++)
	{if(j==i)
		{for(int k=0;k<4;k++)
		{J[i][j]+=-(m[i/2][k]*e[k]-n[i/2][k]*f[k]);}
		J[i][j]+=-m[i/2][i/2]*e[i/2]-n[i/2][i/2]*f[i/2];
				}
        else if(j==i+1)
		{  for(int k=0;k<4;k++)
		{ J[i][j]+=-(m[i/2][k]*f[k]+n[i/2][k]*e[k]);}
		J[i][j]+=n[i/2][i/2]*e[i/2]-m[i/2][i/2]*f[i/2];}
		else if (j%2==0)
		{J[i][j]+=-(m[i/2][j/2]*e[i/2]+n[i/2][j/2]*f[i/2]);}
		else 
		{  J[i][j]+=n[i/2][j/2]*e[i/2]-m[i/2][j/2]*f[i/2];}}}
		else if(i%2!=0&&i<4)
		{
			for(int j=0;j<6;j++)
			{if(j==i)
			{for(int k=0;k<4;k++)
			{J[i][j]+=-(m[i/2][k]*e[k]-n[i/2][k]*f[k]);}
			J[i][j]+=m[i/2][i/2]*e[i/2]+n[i/2][i/2]*f[i/2];}
		else if(j==i-1)
			{for(int k=0;k<4;k++)
			{J[i][j]+=(m[i/2][k]*f[k]+n[i/2][k]*e[k]);}
			J[i][j]+=n[i/2][i/2]*e[i/2]-m[i/2][i/2]*f[i/2];}
		else if(j%2==0)
			{J[i][j]+=n[i/2][j/2]*e[i/2]-m[i/2][j/2]*f[i/2];}
		else 
			{J[i][j]+=m[i/2][j/2]*e[i/2]+n[i/2][j/2]*f[i/2];}} }
		else {for(int j=0;j<6;j++)
		{ if(j==i)
		{      J[i][j]+=-2*f[i/2];}
		else if(j==i-1)
		{      J[i][j]+=-2*e[i/2];}
		else 
		{ J[i][j]=0;     }}}
	}
cout<<"雅可比行列式"<<endl;
for(i=0;i<=5;i++)
{for(j=0;j<=5;j++)
	   cout<<J[i][j]<<"     ";}
double S[6][7];
for(i=0;i<=5;i++)
   for(j=0;j<=5;j++)
	  S[i][j]=J[i][j];
 S[0][6]=p[0];
 S[1][6]=q[0];
 S[2][6]=p[1];
 S[3][6]=q[1];
 S[4][6]=p[2];
 S[5][6]=vv;
 cout<<endl;
cout<<"增廣矩陣"<<endl;
for(i=0;i<=5;i++)
{ for(j=0;j<=6;j++)
cout<<S[i][j]<<"   ";
cout<<endl;}
int l=0,h,z;
double t,var,
      v3[6];
while(l!=5)
{if(S[l][l]==0)
for(h=l+1;h<6;h++)
{if(S[h][l]!=0)
for( z=l;z<7;z++)
{   var=S[l][z];
	S[l][z]=S[h][z];
	S[h][z]=var;} 
break;}
for(i=l+1;i<6;i++)
{ 
	t=S[i][l];
	for(j=l;j<7;j++)
	{S[i][j]+=(-t/S[l][l]*S[l][j]);}
}
	l++;
}

l=5;
while(l!=0)
{	for(i=l-1;i>=0;i--)
{ t=S[i][l];
S[i][6]+=(-t/S[l][l]*S[l][6]);}
l--;}
for(h=0;h<6;h++)
{v3[h]=-S[h][6]/S[h][h];}
cout<<"修正方程的解"<<endl;
for(h=0;h<=5;h++)
cout<<v3[h]<<"  ";
e[0]=e[0]+v3[0];
e[1]=e[1]+v3[2];
e[2]=e[2]+v3[4];
f[0]=f[0]+v3[1];
f[1]=f[1]+v3[3];
f[2]=f[2]+v3[5];
cout<<endl;
cout<<"節(jié)點電壓"<<endl;
cout<<e[0]<<"-j"<<-f[0]<<' '<<e[1]<<"-j"<<-f[1]<<' '<<e[2]<<"+j"<<f[2]<<endl;}

cout<<" 節(jié)點電壓結(jié)果"<<endl;
cout<<"V1="<<e[0]<<"-j"<<-f[0]<<"  "<<"V2="<<e[1]<<"-j"<<-f[1]<<"  "<<"V3="<<e[2]<<"+j"<<f[2]<<endl;
}
}