1、直角坐標(biāo)系下牛頓法潮流計(jì)算1 計(jì)算原理1.1 節(jié)點(diǎn)導(dǎo)納矩陣的形成在圖1（a）的簡(jiǎn)單電力系統(tǒng)中，若略去變壓器的勵(lì)磁功率和線路電容，負(fù)荷用阻抗表示，便可以得到一個(gè)有5個(gè)節(jié)點(diǎn)（包括零電位點(diǎn)）和7條支路的等值網(wǎng)絡(luò)，如圖1（b）所示。將接于節(jié)點(diǎn)1和4的電勢(shì)源和阻抗的串聯(lián)組合變換成等值的電流源和導(dǎo)納的并聯(lián)組合，變得到圖（c）的等值網(wǎng)絡(luò)，其中和分別稱為節(jié)點(diǎn)1和4的注入電流源。 圖1 電力系統(tǒng)及其網(wǎng)絡(luò)以零電位點(diǎn)作為計(jì)算節(jié)點(diǎn)電壓的參考點(diǎn)，根據(jù)基爾霍夫定律，可以寫出4個(gè)獨(dú)立節(jié)點(diǎn)的電流平衡方程如下 （1）上述方程組經(jīng)過(guò)整理可以寫成 （2）式中，；。一般的，對(duì)于有個(gè)獨(dú)立節(jié)點(diǎn)的網(wǎng)絡(luò)，可以列寫個(gè)節(jié)點(diǎn)方程 （3）也可以用矩

2、陣寫成 （4）或縮寫為 （5）矩陣稱為節(jié)點(diǎn)導(dǎo)納矩陣。它的對(duì)角線元素稱為節(jié)點(diǎn)的自導(dǎo)納，其值等于接于節(jié)點(diǎn)的所有支路導(dǎo)納之和。非對(duì)角線元素稱為節(jié)點(diǎn)、 間的互導(dǎo)納，它等于直接接于節(jié)點(diǎn)、間的支路導(dǎo)納的負(fù)值。若節(jié)點(diǎn)、間不存在直接支路，則有。由此可知節(jié)點(diǎn)導(dǎo)納矩陣是一個(gè)稀疏的對(duì)稱矩陣。1.2 直角坐標(biāo)系下牛頓法潮流計(jì)算 采用直角坐標(biāo)時(shí)，節(jié)點(diǎn)電壓可表示為導(dǎo)納矩陣元素則表示為將上述表示式代入的右端，展開并分出實(shí)部和虛部，便得 （6）假定系統(tǒng)中的第1，2，3，m號(hào)節(jié)點(diǎn)為PQ節(jié)點(diǎn)，第i個(gè)節(jié)點(diǎn)的給定功率設(shè)為和，對(duì)對(duì)該節(jié)點(diǎn)可列寫方程 （i=1,2，m） （7）假定系統(tǒng)中的第m+1，m+2，n-1號(hào)節(jié)點(diǎn)為PV節(jié)點(diǎn)，則對(duì)其中

3、每一個(gè)節(jié)點(diǎn)可以列寫方程 （i=m+1，m+2，n-1） （8）第n號(hào)節(jié)點(diǎn)為平衡點(diǎn)，其電壓是給定的，故不參加迭代。式（7）和式（8）總共包含了2（n-1）個(gè)方程，待求的變量有也是2（n-1）個(gè)。我們還可看到，方程（7）和式（8）已經(jīng)具備了方程組的形式。因此，不難寫出如下的修正方程式 （9）式中上述方程中雅克比矩陣的各元素，可以對(duì)式（7）和式（8）求偏導(dǎo)數(shù)獲得。當(dāng)時(shí) （10）當(dāng)時(shí) （11）修正方程式（11-48）還可以寫成分塊矩陣的形式 （12）式中，和都是二維列向量；是介方陣。對(duì)于PQ節(jié)點(diǎn) （13）對(duì)于PV節(jié)點(diǎn) （14）從表達(dá)式（1-7）（1-11）可以看到，雅克比矩陣有以下特點(diǎn)：（1） 雅克比

4、矩陣各元素都是節(jié)點(diǎn)電壓的函數(shù)，它們的數(shù)值將在迭代過(guò)程中不斷的改變。（2） 雅克比矩陣的子塊中的元素的表達(dá)式只用到導(dǎo)納矩陣中的對(duì)應(yīng)元素。若，則必有。因此，式（1-9）式中分塊形式的雅克比矩陣同節(jié)點(diǎn)導(dǎo)納矩陣一樣稀疏，修正方程的求解同樣可以用稀疏矩陣的求解技巧。（3） 無(wú)論在式（1-6）或式（1-9）中雅克比矩陣的元素或子塊都不具有對(duì)稱性。用牛頓-拉夫遜法計(jì)算潮流的流程：首先要輸入網(wǎng)絡(luò)的原始數(shù)據(jù)以及各節(jié)點(diǎn)的給定值并形成節(jié)點(diǎn)導(dǎo)納矩陣。輸入節(jié)點(diǎn)電壓初值和，置迭代計(jì)數(shù)k=0。然后開始進(jìn)入牛頓法的迭代過(guò)程。在進(jìn)行第k+1次迭代時(shí)，其計(jì)算步驟如下：（1） 按上一次迭代計(jì)算出的節(jié)點(diǎn)電壓值和，利用式（7）和式（8

5、）計(jì)算各類節(jié)點(diǎn)的不平衡量、和。（2） 按條件校驗(yàn)收斂，即 如果收斂，迭代到此結(jié)束，轉(zhuǎn)入計(jì)算各線路潮流和平衡節(jié)點(diǎn)的功率，并打印輸出計(jì)算結(jié)果。不收斂則繼續(xù)計(jì)算。（3）利用式（10）和式（11）計(jì)算雅克比矩陣的各元素。（4）解修正方程式（7）求節(jié)點(diǎn)電壓的修正量和。（5）修正各節(jié)點(diǎn)的電壓（6）迭代計(jì)數(shù)加1，返回第一步繼續(xù)迭代過(guò)程。2 計(jì)算過(guò)程21 計(jì)算程序框圖圖2 潮流計(jì)算程序框圖2.2 節(jié)點(diǎn)導(dǎo)納矩陣的形成圖3 節(jié)點(diǎn)導(dǎo)納矩陣根據(jù)第一章節(jié)點(diǎn)導(dǎo)納矩陣的形成原理，設(shè)計(jì)題目可化為如圖含有4節(jié)點(diǎn)的形式。采用直角坐標(biāo)時(shí)，節(jié)點(diǎn)電壓可表示為導(dǎo)納矩陣元素則表示為/*計(jì)算導(dǎo)納矩陣*G11=z12r/(z12r*z12r+

6、z12m*z12m)+k*k*z13r/(z13r*z13r+z13m*z13m)+z14r/(z14r*z14r+z14m*z14m);B11=-z12m/(z12r*z12r+z12m*z12m)-k*k*z13m/(z13r*z13r+z13m*z13m)-z14m/(z14r*z14r+z14m*z14m)+y140+y120;G22=z12r/(z12r*z12r+z12m*z12m)+z24r/(z24r*z24r+z24m*z24m);B22=-z12m/(z12r*z12r+z12m*z12m)-z24m/(z24r*z24r+z24m*z24m)+y240+y120;G33=

7、z13r/(z13r*z13r+z13m*z13m);B33=-z13m/(z13r*z13r+z13m*z13m);G44=z14r/(z14r*z14r+z14m*z14m)+z24r/(z24r*z24r+z24m*z24m);B44=-z14m/(z14r*z14r+z14m*z14m)-z24m/(z24r*z24r+z24m*z24m)+y240+y140;G12=G21=-z12r/(z12r*z12r+z12m*z12m);B12=B21=z12m/(z12r*z12r+z12m*z12m);G13=G31=-k*z13r/(z13r*z13r+z13m*z13m);B13=B

8、31=k*z13m/(z13r*z13r+z13m*z13m);G14=G41=-z14r/(z14r*z14r+z14m*z14m);B14=B41=z14m/(z14r*z14r+z14m*z14m);G23=G32=0.0;B23=B32=0.0;G24=G42=-z24r/(z24r*z24r+z24m*z24m);B24=B42=z24m/(z24r*z24r+z24m*z24m);G34=G43=0.0;B34=B43=0.0;for(i=1;i5;i+)for(j=1;j5;j+)printf(%f+%fj,Gij,Bij);printf( );printf(n);/形成節(jié)點(diǎn)導(dǎo)納

9、矩陣將上述表示式代入的右端，展開并分出實(shí)部和虛部，便得2.3 計(jì)算各節(jié)點(diǎn)不平衡量假定系統(tǒng)中的第1，2，3，m號(hào)節(jié)點(diǎn)為PQ節(jié)點(diǎn)，第i個(gè)節(jié)點(diǎn)的給定功率設(shè)為和，對(duì)對(duì)該節(jié)點(diǎn)可列寫方程 （i=1,2，m） 假定系統(tǒng)中的第m+1，m+2，n-1號(hào)節(jié)點(diǎn)為PV節(jié)點(diǎn)，則對(duì)其中每一個(gè)節(jié)點(diǎn)可以列寫方程 （i=m+1，m+2，n-1） 第n號(hào)節(jié)點(diǎn)為平衡點(diǎn)，其電壓是給定的，故不參加迭代，其計(jì)算程序如下：/計(jì)算各節(jié)點(diǎn)不平衡量loop1:printf(迭代次數(shù)k1=%dn,k1);for (i=1;i5;i+)float a=0,b=0;for(j=1;j5;j+)a+=Gij*ej-Bij*fj; b+=Gij*fj+B

10、ij*ej;Pi=Psi-(ei*a+fi*b);/計(jì)算有功功率的增量Qi=Qsi-(fi*a-ei*b);/計(jì)算無(wú)功功率的增量V32=V3s*V3s-e3*e3;printf(有功功率增量P1=%f,P1); printf( ,);printf(有功功率增量P2=%f,P2); printf( ,);printf(有功功率增量P3=%f,P3);printf(無(wú)功功率增量Q1=%f,Q1); printf( ,);printf(無(wú)功功率增量Q2=%f,Q2); printf( ,);printf(電壓增量V32=%f,V32);printf(n);2.4 雅克比矩陣計(jì)算上述方程中雅克比矩陣的

11、各元素，可以對(duì)計(jì)算各點(diǎn)不平衡量得公式中求偏導(dǎo)數(shù)獲得。當(dāng)時(shí)當(dāng)時(shí)以下為程序，只列出j=1的情況/*形成雅克比矩陣*for(j=1;j4;j+)if(1=j)float c=0,d=0;int m;for(m=1;m5;m+)c+=G1m*em-B1m*fm; d+=G1m*fm+B1m*em;J1*N-1j*N-1=-c-G1j*e1-B1j*f1;J1*N-1j*N=-d+B1j*e1-G1j*f1;J1*Nj*N-1=d+B1j*e1-G1j*f1;J1*Nj*N=-c+G1j*e1+B1j*f1;elseJ1*N-1j*N-1=-G1j*e1-B1j*f1; J1*Nj*N=G1j*e1-B

12、1j*f1; J1*N-1j*N=B1j*e1-G1j*f1; J1*Nj*N-1=B1j*e1-G1j*f1;2.5 LU分解法求修正方程LU分解,又稱Gauss消去法,可把任意方陣分解成下三角矩陣的基本變換形式(行交換)和上三角矩陣的乘積。其數(shù)學(xué)表達(dá)式為:A=LU。其中L為下三角矩陣的基本變換形式,U為上三角矩陣。若有矩陣Ax=b 把矩陣LU分解，求AX=b的問(wèn)題就等價(jià)于求出A=LU后：因?yàn)長(zhǎng)y=b可求y，再因?yàn)閁x=y，可求出x。原始的求法x=A（-1）*b,某些情況下，如果矩陣A中的數(shù)非常小，我認(rèn)為不是因?yàn)榇髷?shù)除以小數(shù)誤差大么，1/A算出的誤差會(huì)很大。但LU可以把A分解成兩個(gè)都比A大的

13、矩陣的乘積，1/L的誤差比1/A小的多。求修正方程的程序如下：/*計(jì)算修正方程*for(i=1;iM;i+)Lii=1;for(i=1;iM;i+)U1i=J1i;Li1=Ji1/U11;for(n=2;nM;n+)for(j=n;jM;j+)sigma1=0;for(s=0;s=n-1;s+)sigma1+=Lns*Usj;Unj=Jnj-sigma1;for(i=n;iM;i+)sigma2=0;for(s=0;s=n-1;s+)sigma2+=Lis*Usn;Lin=(Jin-sigma2)/Unn;b1=P1;b2=Q1;b3=P2;b4=Q2;b5=P3;b6=V32;for(i=1

14、;iM;i+)sigma1=0;for(n=1;n=1;i-)sigma2=0;for(n=i+1;nM;n+)sigma2+=Uin*xn;xi=(yi-sigma2)/Uii;xe1=-x1;xe2=-x3;xe3=-x5;xf1=-x2;xf2=-x4;xf3=-x6;printf(節(jié)點(diǎn)電壓:n);for(i=1;i4;i+)ei+=xei; fi+=xfi;for(i=1;i4;i+)printf(e%d=,i); printf(%f,ei); printf( ,);for(i=1;i4;i+)printf(f%d=,i); printf(%f,fi); printf( ,);2.6

15、計(jì)算網(wǎng)絡(luò)中功率分布最后要計(jì)算出平衡節(jié)點(diǎn)的功率和網(wǎng)絡(luò)中的功率分布。輸電線路功率的計(jì)算公式如下：3 程序運(yùn)行結(jié)果經(jīng)過(guò)C語(yǔ)言程序運(yùn)行，可得以下結(jié)果。圖4 程序運(yùn)行結(jié)果4 小結(jié)與心得體會(huì)本文編寫的牛頓法潮流計(jì)算可以滿足設(shè)計(jì)要求，計(jì)算出網(wǎng)絡(luò)功率分布，并且能夠有屏幕顯示，極大地方便了人機(jī)交互。這次電力系統(tǒng)分析課程設(shè)計(jì)歷時(shí)一個(gè)星期，在這一星期的日子里，不僅鞏固了以前所學(xué)過(guò)的知識(shí)，而且學(xué)到了很多在書本上所沒(méi)有學(xué)到過(guò)的知識(shí)。在這次的課程設(shè)計(jì)中，對(duì)于C語(yǔ)言的各種功能終于有了一個(gè)比較全面和具體的認(rèn)識(shí)，在親自動(dòng)手編寫程序的過(guò)程中，發(fā)現(xiàn)了很多讀程序時(shí)不能發(fā)現(xiàn)的漏洞。在編程之余，對(duì)于電力系統(tǒng)分析的計(jì)算方法也有了新的了解，

16、諸如潮流計(jì)算，節(jié)點(diǎn)導(dǎo)納法等。并且在設(shè)計(jì)過(guò)程中需要借組同學(xué)們的智慧，以及書籍，網(wǎng)絡(luò)解決書本外的問(wèn)題，在與同學(xué)討論的過(guò)程中收獲是非常大的，比如使用LU分解法求修正方程的方法就是經(jīng)過(guò)討論和查閱資料學(xué)習(xí)并使用的。所以光靠我自己的力量是很難完成任務(wù)的，也明白了三人行必有我?guī)煹牡览?。通過(guò)這次課程設(shè)計(jì)使我更加體會(huì)到了理論與實(shí)際相結(jié)合的重要性，只有理論知識(shí)是遠(yuǎn)遠(yuǎn)不夠的，在實(shí)踐中可能會(huì)遇到各種各樣的問(wèn)題，不多經(jīng)歷就無(wú)法感受到這一點(diǎn)。要在實(shí)踐中提高自己的動(dòng)手能力和解決問(wèn)題的能力，從而學(xué)以致用。最后，要感謝我的老師們的指導(dǎo)，他們不僅教會(huì)我專業(yè)必須掌握的知識(shí)技能，而且也使我懂得自主學(xué)習(xí)，持之以恒的道理，為將來(lái)的學(xué)習(xí)和

17、工作夯實(shí)基礎(chǔ)。參考文獻(xiàn)1 何仰贊等.電力系統(tǒng)分析上冊(cè)M武漢：華中理工大學(xué)出版社.2 何仰贊等.電力系統(tǒng)分析下冊(cè)M武漢：華中理工大學(xué)出版社.3 諸俊偉等.電力系統(tǒng)分析M.北京：中國(guó)電力出版社，1995.4 周全仁等.電網(wǎng)計(jì)算與程序設(shè)計(jì)M.長(zhǎng)沙：湖南科學(xué)技術(shù)出版社,1983.5丁化成.單片機(jī)應(yīng)用技術(shù)A.北京:北京航空航天大學(xué)出版社,2000.附錄C語(yǔ)言程序：#include #include #include #define N 2#define M 7main()float z12r,z12m,y120,z13r,z13m,k,z14r,z14m,y140,z24r,z24m,y240,G55,

18、B55,J77;float e5=0,1,1,1.1,1.05,f5=0,P5,Q5,Ps5=0,-0.30,-0.55,0.5,xe4,xf4;float Qs5=0,-0.18,-0.13,V3s=1.1,V4S=1.05;float V32,max,P4,Q4;float a1=0,b1=0;int i,j,n,s,k1=0;float LMM=0,UMM=0,sigma1,sigma2,bM,yM,xM;float p12,p13,p14,p21,p24,p31,p41,p42,q12,q13,q14,q21,q24,q31,q41,q42;printf(請(qǐng)輸入z12的實(shí)部和虛部n);

19、 /輸入電路中的阻抗scanf(%f%f,&z12r,&z12m);printf(請(qǐng)輸入z13的實(shí)部和虛部n);scanf(%f%f,&z13r,&z13m);printf(請(qǐng)輸入z14的實(shí)部和虛部n);scanf(%f%f,&z14r,&z14m);printf(請(qǐng)輸入z24的實(shí)部和虛部n);scanf(%f%f,&z24r,&z24m);printf(請(qǐng)輸入y120的值n);scanf(%f,&y120);printf(請(qǐng)輸入y140的值n);scanf(%f,&y140);printf(請(qǐng)輸入y240的值n);scanf(%f,&y240);printf(請(qǐng)輸入變比k的值n);scanf

20、(%f,&k);/*計(jì)算導(dǎo)納矩陣*G11=z12r/(z12r*z12r+z12m*z12m)+k*k*z13r/(z13r*z13r+z13m*z13m)+z14r/(z14r*z14r+z14m*z14m);B11=-z12m/(z12r*z12r+z12m*z12m)-k*k*z13m/(z13r*z13r+z13m*z13m)-z14m/(z14r*z14r+z14m*z14m)+y140+y120;G22=z12r/(z12r*z12r+z12m*z12m)+z24r/(z24r*z24r+z24m*z24m);B22=-z12m/(z12r*z12r+z12m*z12m)-z24m

21、/(z24r*z24r+z24m*z24m)+y240+y120;G33=z13r/(z13r*z13r+z13m*z13m);B33=-z13m/(z13r*z13r+z13m*z13m);G44=z14r/(z14r*z14r+z14m*z14m)+z24r/(z24r*z24r+z24m*z24m);B44=-z14m/(z14r*z14r+z14m*z14m)-z24m/(z24r*z24r+z24m*z24m)+y240+y140;G12=G21=-z12r/(z12r*z12r+z12m*z12m);B12=B21=z12m/(z12r*z12r+z12m*z12m);G13=G3

22、1=-k*z13r/(z13r*z13r+z13m*z13m);B13=B31=k*z13m/(z13r*z13r+z13m*z13m);G14=G41=-z14r/(z14r*z14r+z14m*z14m);B14=B41=z14m/(z14r*z14r+z14m*z14m);G23=G32=0.0;B23=B32=0.0;G24=G42=-z24r/(z24r*z24r+z24m*z24m);B24=B42=z24m/(z24r*z24r+z24m*z24m);G34=G43=0.0;B34=B43=0.0;for(i=1;i5;i+)for(j=1;j5;j+)printf(%f+%fj

23、,Gij,Bij);printf( );printf(n);/形成節(jié)點(diǎn)導(dǎo)納矩陣/*printf(n);/*/計(jì)算各節(jié)點(diǎn)不平衡量loop1:printf(迭代次數(shù)k1=%dn,k1);for (i=1;i5;i+)float a=0,b=0;for(j=1;jfabs(P2)?fabs(P1):fabs(P2);max=maxfabs(P3)?max:fabs(P3);max=maxfabs(Q1)?max:fabs(Q1);max=maxfabs(Q2)?max:fabs(Q2);max=maxfabs(V32)?max:fabs(V32);printf(max=%fn,max);/*whil

24、e (max0.00001)/*形成雅克比矩陣*for(j=1;j4;j+)if(1=j)float c=0,d=0;int m;for(m=1;m5;m+)c+=G1m*em-B1m*fm; d+=G1m*fm+B1m*em;J1*N-1j*N-1=-c-G1j*e1-B1j*f1;J1*N-1j*N=-d+B1j*e1-G1j*f1;J1*Nj*N-1=d+B1j*e1-G1j*f1;J1*Nj*N=-c+G1j*e1+B1j*f1;elseJ1*N-1j*N-1=-G1j*e1-B1j*f1; J1*Nj*N=G1j*e1-B1j*f1; J1*N-1j*N=B1j*e1-G1j*f1;

25、 J1*Nj*N-1=B1j*e1-G1j*f1;for(j=1;j4;j+)if(2=j)float c=0,d=0;int m;for(m=1;m5;m+)c+=G2m*em-B2m*fm; d+=G2m*fm+B2m*em;J2*N-1j*N-1=-c-G2j*e2-B2j*f2;J2*N-1j*N=-d+B2j*e2-G2j*f2;J2*Nj*N-1=d+B2j*e2-G2j*f2;J2*Nj*N=-c+G2j*e2+B2j*f2;elseJ2*N-1j*N-1=-G2j*e2-B2j*f2;J2*Nj*N=G2j*e2-B2j*f2;J2*N-1j*N=B2j*e2-G2j*f2;J

26、2*Nj*N-1=B2j*e2-G2j*f2;for(j=1;j4;j+)if(3=j)float c=0,d=0;int m;for(m=1;m5;m+)c+=G3m*em-B3m*fm; d+=G3m*fm+B3m*em;J3*N-1j*N-1=-c-G3j*e3-B3j*f3;J3*N-1j*N=-d+B3j*e3-G3j*f3;J3*Nj*N-1=-2*e3;J3*Nj*N=-2*f3;elseJ3*N-1j*N-1=-G3j*e3-B3j*f3; J3*N-1j*N=B3j*e3-G3j*f3; J3*Nj*N-1=0; J3*Nj*N=0;printf(雅克比矩陣是:n);for(

27、i=1;i7;i+)for(j=1;j7;j+)printf(%f,Jij);printf( );printf(n);/*計(jì)算修正方程*for(i=1;iM;i+)Lii=1;for(i=1;iM;i+)U1i=J1i;Li1=Ji1/U11;for(n=2;nM;n+)for(j=n;jM;j+)sigma1=0;for(s=0;s=n-1;s+)sigma1+=Lns*Usj;Unj=Jnj-sigma1;for(i=n;iM;i+)sigma2=0;for(s=0;s=n-1;s+)sigma2+=Lis*Usn;Lin=(Jin-sigma2)/Unn;b1=P1;b2=Q1;b3=P

28、2;b4=Q2;b5=P3;b6=V32;for(i=1;iM;i+)sigma1=0;for(n=1;n=1;i-)sigma2=0;for(n=i+1;nM;n+)sigma2+=Uin*xn;xi=(yi-sigma2)/Uii;xe1=-x1;xe2=-x3;xe3=-x5;xf1=-x2;xf2=-x4;xf3=-x6;printf(節(jié)點(diǎn)電壓:n);for(i=1;i4;i+)ei+=xei; fi+=xfi;for(i=1;i4;i+)printf(e%d=,i); printf(%f,ei); printf( ,);for(i=1;i4;i+)printf(f%d=,i); pr

29、intf(%f,fi); printf( ,);printf(n);k1=k1+1;goto loop1;for(j=1;j5;j+)a1+=G4j*ej-B4j*fj; b1+=G4j*fj+B4j*ej;P4=e4*a1+f4*b1;Q4=f4*a1-e4*b1;printf(P4+Q4=%f+j%f,P4,Q4);printf(n);p12=-2*e1*f1*y120-(e1*(e1-e2)-f1*(f1-f2)*G12+(e1*(f1-f2)+f1*(e1-e2)*B12;q12=-(e1*e1-f1*f1)*y120+(e1*(e1-e2)+f1*(f1-f2)*B12+(e1*(f

30、1-f2)+f1*(e1-e2)*G12;p13=(e1*e1-f1*f1)*k*(k-1)*z13r/(z13r*z13r+z13m*z13m)+2*e1*f1*k*(k-1)*z13m/(z13r*z13r+z13m*z13m)-(e1*(e1-e3)-f1*(f1-f3)*G13+(e1*(f1-f3)+f1*(e1-e3)*B13;q13=(e1*e1-f1*f1)*k*(k-1)*z13m/(z13r*z13r+z13m*z13m)-2*e1*f1*k*(k-1)*z13r/(z13r*z13r+z13m*z13m)+(e1*(e1-e3)-f1*(f1-f3)*B13-(e1*(f1-f3)+f1*(e1-e3)*G13;p14=-2*e1*f1*y140-(e1*(e1-e4)-f1*(f1-f4)*G14+(e1*(f1-f4)+f1*(e1-e4)*B14;q14=-(e1*e1-f1*f1)*y140+(e1*(e1-e4)-f1*(f1-f4)*B14+(e1*(f1-f4)+f1*(e1-e4)*G14;p21=-2*e2*f2*y120-(e2*(e2-e1)-f2*(f2-f1)*G21+(e2*(f2-f1)+f2*(e2-e1)*B21;q21=-(e2*