adfunc.cpp 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835
  1. /********* Functions needed for the main program *********/
  2. #include "stdafx.h"
  3. extern int Nc, Die, Rec, Cir, Tri, Ell, Pol, Gnd, Tra;
  4. extern Matrix Gnddi, Recdi, Cirdi, Tradix, Tradiy;
  5. extern Vector ag, bg, cg;
  6. extern Vector a, b, c;
  7. extern Vector at, bt, ct;
  8. extern Vector Circum;
  9. extern Vector Diec, Recc, Circ, Tric, Ellc, Polc, Gndc, Trac;
  10. extern double freq, sqrtoms, sigma, omega, oms, omegamu;
  11. extern int matr, Mt, Np;
  12. extern int Ns, Nsr, Ngx, Ngy, Nxt, Nyt, Ngr, Nrt;
  13. extern int Nh, Nc4, Nit, J, Nw, Nwx, Nwy, Nws;
  14. extern double hsrw, EPS;
  15. extern double Pi, Pih, Pit, wn, mu, cf;
  16. extern double step_w, hw, hwh, mmax;
  17. extern Complex Ic;
  18. /*****************************
  19. double sqv( double x ) {
  20. return( x*x );
  21. }
  22. double csft( double x ) {
  23. if( x < 0.0 ) {
  24. return( x + 1.0 );
  25. } else if( x >= 1.0 ) {
  26. return( x - 1.0 );
  27. } else {
  28. return( x );
  29. }
  30. }
  31. ********************/
  32. /* To calculate the length of a side */
  33. /********************
  34. double sidelen(double x1, double y1, double x2, double y2 ) {
  35. return(sqrt(sqv(x1-x2)+sqv(y1-y2)));
  36. }
  37. *************************/
  38. /* To calculate the area of trapezoid with the coodinates */
  39. double traparea( int i ) {
  40. int j, j1, flag;
  41. double side1, side2, height;
  42. double x[4], y[4], xmax;
  43. for(j = 0; j < 4; j ++) {
  44. y[j] = Tradix(i,j);
  45. }
  46. height = fabs(Tradiy(i,0)-Tradiy(i,2));
  47. for(j = 0; j < 4; j ++) {
  48. xmax = -100;
  49. for(j1 = 0; j1 < 4; j1 ++) {
  50. if(y[j1] > xmax) {
  51. xmax = y[j1];
  52. flag = j1;
  53. }
  54. }
  55. x[j] = xmax;
  56. y[flag] = -200;
  57. }
  58. side1 = fabs(x[0]-x[3]);
  59. side2 = fabs(x[1]-x[2]);
  60. return((side1+side2)*height/2.0);
  61. }
  62. /* Determine the coordinates of rectangular shape */
  63. int gndgeom( int i, double& ti, double& xt, double& yt ) {
  64. register double t = ti*Circum[i];
  65. if(t < 0.0 || t > Circum[i]) {
  66. cerr << " t is out of range in gndgeom !!! " << endl;
  67. exit( 1 );
  68. }
  69. if(t >= 0.0 && t <= ag[i]) {
  70. xt = Gnddi(i,0)+t;
  71. yt = Gnddi(i,1);
  72. return( 1 );
  73. } else if(t > ag[i] && t <= bg[i]) {
  74. xt = Gnddi(i,0)+Gnddi(i,2);
  75. yt = t-ag[i]+Gnddi(i,1);
  76. return( 2 );
  77. } else if(t > bg[i] && t <= cg[i]) {
  78. xt = Gnddi(i,2)-t+bg[i]+Gnddi(i,0);
  79. yt = Gnddi(i,3)+Gnddi(i,1);
  80. return( 3 );
  81. } else if(t > cg[i] && t <= cg[i]+Gnddi(i,3)) {
  82. xt = Gnddi(i,0);
  83. yt = Gnddi(i,3)-t+cg[i]+Gnddi(i,1);
  84. return( 4 );
  85. }
  86. return( 0 );
  87. }
  88. int recgeom( int i, double& ti, double& xt, double& yt ) {
  89. register double t = ti*Circum[i];
  90. int ii = i - Gnd;
  91. if(t < 0.0 || t > Circum[i]) {
  92. cerr << " t is out of range in recgeom !!! " << endl;
  93. exit( 1 );
  94. }
  95. if(t >= 0.0 && t <= a[ii]) {
  96. xt = Recdi(ii,0)+t;
  97. yt = Recdi(ii,1);
  98. return( 1 );
  99. } else if(t > a[ii] && t <= b[ii]) {
  100. xt = Recdi(ii,0)+Recdi(ii,2);
  101. yt = t-a[ii]+Recdi(ii,1);
  102. return( 2 );
  103. } else if(t > b[ii] && t <= c[ii]) {
  104. xt = Recdi(ii,2)-t+b[ii]+Recdi(ii,0);
  105. yt = Recdi(ii,3)+Recdi(ii,1);
  106. return( 3 );
  107. } else if(t > c[ii] && t <= c[ii]+Recdi(ii,3)) {
  108. xt = Recdi(ii,0);
  109. yt = Recdi(ii,3)-t+c[ii]+Recdi(ii,1);
  110. return( 4 );
  111. }
  112. return( 0 );
  113. }
  114. int cirgeom(int i, double& ti, double& xt, double& yt, double& al) {
  115. register double t = ti*Circum[i];
  116. int ii = i-Gnd-Rec;
  117. if( t < 0.0 || t > Circum[i] ) {
  118. cerr << " t is out of range in cirgeom !!! " << endl;
  119. exit( 1 );
  120. } else {
  121. al = t / Cirdi(ii,2);
  122. xt = Cirdi(ii,0) + Cirdi(ii,2) * cos(al);
  123. yt = Cirdi(ii,1) + Cirdi(ii,2) * sin(al);
  124. return( 5 );
  125. }
  126. return( 0 );
  127. }
  128. int trageom( int i, double& ti, double& xt, double& yt ) {
  129. register double t = ti*Circum[i];
  130. double tside;
  131. int ii=i-Gnd-Rec-Cir;
  132. if( t < 0.0 || t > Circum[i] ) {
  133. cerr << " t is out of range in trageom !!! " << endl;
  134. exit( 1 );
  135. }
  136. if( t >= 0.0 && t < at[ii] ) {
  137. tside = sidelen(Tradix(ii,0),Tradiy(ii,0),Tradix(ii,1),Tradiy(ii,1));
  138. xt = Tradix(ii,0)+t*(Tradix(ii,1)-Tradix(ii,0))/tside;
  139. yt = Tradiy(ii,0)+t*(Tradiy(ii,1)-Tradiy(ii,0))/tside;
  140. return( 6 );
  141. } else if( t > at[ii] && t <= bt[ii] ) {
  142. t -= at[ii];
  143. tside = sidelen(Tradix(ii,1),Tradiy(ii,1),Tradix(ii,2),Tradiy(ii,2));
  144. xt = Tradix(ii,1)+t*(Tradix(ii,2)-Tradix(ii,1))/tside;
  145. yt = Tradiy(ii,1)+t*(Tradiy(ii,2)-Tradiy(ii,1))/tside;
  146. return( 7 );
  147. } else if( t > bt[ii] && t <= ct[ii] ) {
  148. t -= bt[ii];
  149. tside = sidelen(Tradix(ii,2),Tradiy(ii,2),Tradix(ii,3),Tradiy(ii,3));
  150. xt = Tradix(ii,2)+t*(Tradix(ii,3)-Tradix(ii,2))/tside;
  151. yt = Tradiy(ii,2)+t*(Tradiy(ii,3)-Tradiy(ii,2))/tside;
  152. return( 8 );
  153. } else if( t > ct[ii] && t <= Circum[i]) {
  154. t -= ct[ii];
  155. tside = sidelen(Tradix(ii,3),Tradiy(ii,3),Tradix(ii,0),Tradiy(ii,0));
  156. xt = Tradix(ii,3)+t*(Tradix(ii,0)-Tradix(ii,3))/tside;
  157. yt = Tradiy(ii,3)+t*(Tradiy(ii,0)-Tradiy(ii,3))/tside;
  158. return( 9 );
  159. }
  160. return( 0 );
  161. }
  162. void getsbs( Matrix& uns ) {
  163. for( int i = 0; i < Nw; i++ ) {
  164. uns(0,i) = i;
  165. uns(1,i) = double(i-4.0)*step_w;
  166. uns(2,i) = double(i+7.0)*step_w;
  167. }
  168. }
  169. // determine the coefficients in the gauss quadrature.
  170. void intpnt( double xcnt, Vector& tx, Vector& wx, Vector& ptx, Vector& whx,
  171. double ycnt, Vector& ty, Vector& wy, Vector& pty, Vector& why,
  172. double h ) {
  173. const int nx = tx.dim();
  174. const int ny = ty.dim();
  175. int i;
  176. for( i = 0; i < nx; i++ ){
  177. tx[i] = h*ptx[i]+xcnt;
  178. wx[i] = whx[i]*h;
  179. }
  180. for( i = 0; i < ny; i++ ) {
  181. ty[i] = h*pty[i]+ycnt;
  182. wy[i] = why[i]*h;
  183. }
  184. }
  185. // calculate the G0
  186. Complex kernel_v0( double& t, double& tp, int i, int j ) {
  187. // register Complex res;
  188. double tt, ttp;
  189. tt = csft( t );
  190. ttp = csft( tp );
  191. // res = g0mn( i, j, tt, ttp );
  192. return( g0mn( i, j, tt, ttp ) );
  193. }
  194. // calculate the (partial G0/partial n')
  195. Complex kernel_u0( double& t, double& tp, int i, int j ) {
  196. // register Complex res;
  197. double tt, ttp;
  198. tt = csft( t );
  199. ttp = csft( tp );
  200. // res = g0pmn( i, j, tt, ttp );
  201. return( g0pmn( i, j, tt, ttp ) );
  202. }
  203. // calculate the Gi
  204. Complex kernel_vi( double& t, double& tp, int i, int j ) {
  205. // register Complex res;
  206. double tt, ttp;
  207. tt = csft( t );
  208. ttp = csft( tp );
  209. // res = gimn( i, j, tt, ttp );
  210. return( gimn( i, j, tt, ttp ) );
  211. }
  212. // calculate the (partial Gi/partial n')
  213. Complex kernel_ui( double& t, double& tp, int i, int j ) {
  214. // register Complex res;
  215. double tt, ttp;
  216. tt = csft( t );
  217. ttp = csft( tp );
  218. // res = gipmn( i, j, tt, ttp );
  219. return( gipmn( i, j, tt, ttp ) );
  220. }
  221. // Pi = 4.0*atan(1.0);
  222. static double s_2Pi = 8.0 * atan(1.0);
  223. // calculate G0
  224. double g0(double x, double y, double xp, double yp) {
  225. // double rou2;
  226. // double g0;
  227. double tmp1 = x - xp;
  228. double tmp2 = y - yp;
  229. // rou2 = tmp1 * tmp1 + tmp2 * tmp2;
  230. // g0 = -1.0/(4.0*Pi)*log(tmp1 * tmp1 + tmp2 * tmp2);
  231. return(-1.0/(4.0*Pi)*log(tmp1 * tmp1 + tmp2 * tmp2));
  232. }
  233. // calculate derivative of G0 in term of x
  234. double g0x(double x, double y, double xp, double yp) {
  235. // double rou2;
  236. // double g0x;
  237. double tmp1 = x - xp;
  238. double tmp2 = y - yp;
  239. // rou2 = tmp1 * tmp1 + tmp2 * tmp2;
  240. // g0x = -tmp1/(s_2Pi*(tmp1 * tmp1 + tmp2 * tmp2));
  241. return(-tmp1/(s_2Pi*(tmp1 * tmp1 + tmp2 * tmp2)));
  242. }
  243. // calculate derivative of G0 in term of y
  244. double g0y(double x, double y, double xp, double yp) {
  245. // double rou2;
  246. // double g0y;
  247. double tmp1 = x - xp;
  248. double tmp2 = y - yp;
  249. // rou2 = tmp1 * tmp1 + tmp2 * tmp2;
  250. // g0y = -tmp2/(s_2Pi*(tmp1 * tmp1 + tmp2 * tmp2));
  251. return(-tmp2/(s_2Pi*(tmp1 * tmp1 + tmp2 * tmp2)));
  252. }
  253. // calculate value of G0
  254. double g0mn(int i, int j, double r, double rp) {
  255. double x,y,xp,yp,al,alp;
  256. double g0mn;
  257. if( i < Gnd ) {
  258. gndgeom(i,r,x,y);
  259. } else if( i >= Gnd && i < (Gnd + Rec) ) {
  260. recgeom(i,r,x,y);
  261. } else if( i >= (Gnd + Rec) && i < (Gnd + Rec + Cir) ) {
  262. cirgeom(i,r,x,y,al);
  263. } else if( i >= (Gnd + Rec + Cir) && i < (Gnd + Rec + Cir + Tra) ) {
  264. trageom(i,r,x,y);
  265. }
  266. if( j < Gnd ) {
  267. gndgeom(j,rp,xp,yp);
  268. } else if( j >= Gnd && j < (Gnd + Rec) ) {
  269. recgeom(j,rp,xp,yp);
  270. } else if( j >= (Gnd + Rec) && j < (Gnd + Rec + Cir) ) {
  271. cirgeom(j,rp,xp,yp,alp);
  272. } else if( j >= (Gnd + Rec + Cir) && j < (Gnd + Rec + Cir + Tra) ) {
  273. trageom(j,rp,xp,yp);
  274. }
  275. g0mn = g0(x,y,xp,yp);
  276. return(g0mn);
  277. }
  278. // calculate the derivative
  279. double g0pmn(int i, int j, double r, double rp) {
  280. double x, y, xp, yp, al, alp;
  281. int n, np;
  282. double g0pmn;
  283. Vector tmpv(2);
  284. int ii = j-Gnd-Rec-Cir;
  285. if( i < Gnd ) {
  286. n = gndgeom(i,r,x,y);
  287. } else if( i >= Gnd && i < (Gnd + Rec) ) {
  288. n = recgeom(i,r,x,y);
  289. } else if( i >= (Gnd + Rec) && i < (Gnd + Rec + Cir) ) {
  290. n = cirgeom(i,r,x,y,al);
  291. } else if( i >= (Gnd + Rec + Cir) && i < (Gnd + Rec + Cir + Tra) ) {
  292. n = trageom(i,r,x,y);
  293. }
  294. if( j < Gnd ) {
  295. np = gndgeom(j,rp,xp,yp);
  296. } else if( j >= Gnd && j < (Gnd + Rec) ) {
  297. np = recgeom(j,rp,xp,yp);
  298. } else if( j >= (Gnd + Rec) && j < (Gnd + Rec + Cir) ) {
  299. np = cirgeom(j,rp,xp,yp,alp);
  300. } else if( j >= (Gnd + Rec + Cir) && j < (Gnd + Rec + Cir + Tra) ) {
  301. np = trageom(j,rp,xp,yp);
  302. }
  303. switch( np ) {
  304. case 1:
  305. g0pmn = g0y(x,y,xp,yp);
  306. break;
  307. case 2:
  308. g0pmn = -g0x(x,y,xp,yp);
  309. break;
  310. case 3:
  311. g0pmn = -g0y(x,y,xp,yp);
  312. break;
  313. case 4:
  314. g0pmn = g0x(x,y,xp,yp);
  315. break;
  316. case 5: // circle
  317. g0pmn = -g0x(x,y,xp,yp)*cos(alp) - g0y(x,y,xp,yp)*sin(alp);
  318. break;
  319. case 6:
  320. tmpv[0] = -g0x(x,y,xp,yp);
  321. tmpv[1] = -g0y(x,y,xp,yp);
  322. g0pmn = tmpv * tangle(Tradix(ii,0),Tradiy(ii,0),Tradix(ii,1),Tradiy(ii,1),Tradix(ii,2),Tradiy(ii,2));
  323. break;
  324. case 7:
  325. tmpv[0] = -g0x(x,y,xp,yp);
  326. tmpv[1] = -g0y(x,y,xp,yp);
  327. g0pmn = tmpv * tangle(Tradix(ii,1),Tradiy(ii,1),Tradix(ii,2),Tradiy(ii,2),Tradix(ii,3),Tradiy(ii,3));
  328. break;
  329. case 8:
  330. tmpv[0] = -g0x(x,y,xp,yp);
  331. tmpv[1] = -g0y(x,y,xp,yp);
  332. g0pmn = tmpv * tangle(Tradix(ii,2),Tradiy(ii,2),Tradix(ii,3),Tradiy(ii,3),Tradix(ii,0),Tradiy(ii,0));
  333. break;
  334. case 9: // 6-9 Trapezoid
  335. tmpv[0] = -g0x(x,y,xp,yp);
  336. tmpv[1] = -g0y(x,y,xp,yp);
  337. g0pmn = tmpv * tangle(Tradix(ii,3),Tradiy(ii,3),Tradix(ii,0),Tradiy(ii,0),Tradix(ii,1),Tradiy(ii,1));
  338. break;
  339. default:
  340. cerr << np << endl;
  341. cerr << "unidentified number" << endl;
  342. }
  343. return(g0pmn);
  344. }
  345. // functions about G1(r,rp)
  346. // calculate the real part of Gi
  347. double gir(double x, double y, double xp, double yp) {
  348. extern double ker(double);
  349. double rou;
  350. // double gir;
  351. double tmpx = x - xp;
  352. double tmpy = y - yp;
  353. rou = sqrt(tmpx*tmpx + tmpy*tmpy)*sqrtoms;
  354. // gir = 1.0/(s_2Pi)*ker(rou);
  355. return(1.0/(s_2Pi)*ker(rou));
  356. }
  357. // calculate the derivative of real part of Gi in term of x
  358. double girx(double x, double y, double xp, double yp) {
  359. extern double kerp(double);
  360. double rou; // ,rouoms;
  361. // double girx;
  362. double tmpx = x - xp;
  363. double tmpy = y - yp;
  364. rou = sqrt(tmpx*tmpx + tmpy*tmpy);
  365. // rouoms = sqrtoms*rou;
  366. // girx = sqrtoms*kerp(sqrtoms*rou)*tmpx/(s_2Pi*rou);
  367. return(sqrtoms*kerp(sqrtoms*rou)*tmpx/(s_2Pi*rou));
  368. }
  369. // calculate the derivative of real part of Gi in term of y
  370. double giry(double x, double y, double xp, double yp) {
  371. extern double kerp(double);
  372. double rou; // ,rouoms;
  373. // double giry;
  374. double tmpx = x - xp;
  375. double tmpy = y - yp;
  376. rou = sqrt(tmpx*tmpx + tmpy*tmpy);
  377. // rouoms = rou*sqrtoms;
  378. // giry = sqrtoms*kerp(rou*sqrtoms)*tmpy/(s_2Pi*rou);
  379. return(sqrtoms*kerp(rou*sqrtoms)*tmpy/(s_2Pi*rou));
  380. }
  381. // calculate the value of Gi
  382. double girmn(int i, int j, double r, double rp) {
  383. double x,y,xp,yp,al,alp;
  384. // double girmn;
  385. if( i < Gnd ) {
  386. gndgeom(i,r,x,y);
  387. } else if( i >= Gnd && i < (Gnd + Rec) ) {
  388. recgeom(i,r,x,y);
  389. } else if( i >= (Gnd + Rec) && i < (Gnd + Rec + Cir) ) {
  390. cirgeom(i,r,x,y,al);
  391. } else if( i >= (Gnd + Rec + Cir) && i < (Gnd + Rec + Cir + Tra) ) {
  392. trageom(i,r,x,y);
  393. }
  394. if( j < Gnd ) {
  395. gndgeom(j,rp,xp,yp);
  396. } else if( j >= Gnd && j < (Gnd + Rec) ) {
  397. recgeom(j,rp,xp,yp);
  398. } else if( j >= (Gnd + Rec) && j < (Gnd + Rec + Cir) ) {
  399. cirgeom(j,rp,xp,yp,alp);
  400. } else if( j >= (Gnd + Rec + Cir) && j < (Gnd + Rec + Cir + Tra) ) {
  401. trageom(j,rp,xp,yp);
  402. }
  403. // girmn = gir(x,y,xp,yp);
  404. return(gir(x,y,xp,yp));
  405. }
  406. // calculate the derivative of Gi
  407. double girpmn(int i, int j, double r, double rp) {
  408. double x,y,xp,yp,al,alp;
  409. double girpmn;
  410. int n,np;
  411. Vector tmpv(2);
  412. int ii=j-Gnd-Rec-Cir;
  413. if( i < Gnd ) {
  414. n = gndgeom(i,r,x,y);
  415. } else if( i >= Gnd && i < (Gnd + Rec) ) {
  416. n = recgeom(i,r,x,y);
  417. } else if( i >= (Gnd + Rec) && i < (Gnd + Rec + Cir) ) {
  418. n = cirgeom(i,r,x,y,al);
  419. } else if( i >= (Gnd + Rec + Cir) && i < (Gnd + Rec + Cir + Tra) ) {
  420. n = trageom(i,r,x,y);
  421. }
  422. if( j < Gnd ) {
  423. np = gndgeom(j,rp,xp,yp);
  424. } else if( j >= Gnd && j < (Gnd + Rec) ) {
  425. np = recgeom(j,rp,xp,yp);
  426. } else if( j >= (Gnd + Rec) && j < (Gnd + Rec + Cir) ) {
  427. np = cirgeom(j,rp,xp,yp,alp);
  428. } else if( j >= (Gnd + Rec + Cir) && j < (Gnd + Rec + Cir + Tra) ) {
  429. np = trageom(j,rp,xp,yp);
  430. }
  431. switch( np ) {
  432. case 1:
  433. girpmn = giry(x,y,xp,yp);
  434. break;
  435. case 2:
  436. girpmn = -girx(x,y,xp,yp);
  437. break;
  438. case 3:
  439. girpmn = -giry(x,y,xp,yp);
  440. break;
  441. case 4:
  442. girpmn = girx(x,y,xp,yp);
  443. break;
  444. case 5: // circle
  445. girpmn = -girx(x,y,xp,yp)*cos(alp) - giry(x,y,xp,yp)*sin(alp);
  446. break;
  447. case 6:
  448. tmpv[0] = -girx(x,y,xp,yp);
  449. tmpv[1] = -giry(x,y,xp,yp);
  450. girpmn = tmpv * tangle(Tradix(ii,0),Tradiy(ii,0),Tradix(ii,1),Tradiy(ii,1),Tradix(ii,2),Tradiy(ii,2));
  451. break;
  452. case 7:
  453. tmpv[0] = -girx(x,y,xp,yp);
  454. tmpv[1] = -giry(x,y,xp,yp);
  455. girpmn = tmpv * tangle(Tradix(ii,1),Tradiy(ii,1),Tradix(ii,2),Tradiy(ii,2),Tradix(ii,3),Tradiy(ii,3));
  456. break;
  457. case 8:
  458. tmpv[0] = -girx(x,y,xp,yp);
  459. tmpv[1] = -giry(x,y,xp,yp);
  460. girpmn = tmpv * tangle(Tradix(ii,2),Tradiy(ii,2),Tradix(ii,3),Tradiy(ii,3),Tradix(ii,0),Tradiy(ii,0));
  461. break;
  462. case 9: // 6-9 Trapezoid
  463. tmpv[0] = -girx(x,y,xp,yp);
  464. tmpv[1] = -giry(x,y,xp,yp);
  465. girpmn = tmpv * tangle(Tradix(ii,3),Tradiy(ii,3),Tradix(ii,0),Tradiy(ii,0),Tradix(ii,1),Tradiy(ii,1));
  466. break;
  467. default:
  468. cerr << np << endl;
  469. cerr << "unidentified number" << endl;
  470. }
  471. return(girpmn);
  472. }
  473. // calculate the imaginary part of Gi
  474. double gii(double x, double y, double xp, double yp) {
  475. extern double kei(double);
  476. double rou;
  477. double gii;
  478. rou = sqrt(sqv(x-xp)+sqv(y-yp))*sqrtoms;
  479. gii = -1.0/(s_2Pi)*kei(rou);
  480. return(gii);
  481. }
  482. // calculate the derivative of imaginary part of Gi in term of x
  483. double giix(double x, double y, double xp, double yp) {
  484. extern double keip(double);
  485. double rou,rouoms;
  486. double giix;
  487. rou = sqrt(sqv(x-xp)+sqv(y-yp));
  488. rouoms = rou*sqrtoms;
  489. if(rouoms>0.0)
  490. giix = -sqrtoms*keip(sqrtoms*rou)*(x-xp)/(s_2Pi*rou);
  491. else
  492. giix = 0.0;
  493. return(giix);
  494. }
  495. // calculate the derivative of imaginary part of Gi in term of y
  496. double giiy(double x, double y, double xp, double yp) {
  497. extern double keip(double);
  498. double rou,rouoms;
  499. double giiy;
  500. rou = sqrt(sqv(x-xp)+sqv(y-yp));
  501. rouoms = rou*sqrtoms;
  502. if(rouoms > 0.0)
  503. giiy = -sqrtoms*keip(sqrtoms*rou)*(y-yp)/(s_2Pi*rou);
  504. else
  505. giiy = 0.0;
  506. return(giiy);
  507. }
  508. // calculate the value of imaginary part of Gi
  509. double giimn(int i, int j, double r, double rp) {
  510. double x,y,xp,yp,al,alp;
  511. double giimn;
  512. if( i < Gnd ) {
  513. gndgeom(i,r,x,y);
  514. } else if( i >= Gnd && i < (Gnd + Rec) ) {
  515. recgeom(i,r,x,y);
  516. } else if( i >= (Gnd + Rec) && i < (Gnd + Rec + Cir) ) {
  517. cirgeom(i,r,x,y,al);
  518. } else if( i >= (Gnd + Rec + Cir) && i < (Gnd + Rec + Cir + Tra) ) {
  519. trageom(i,r,x,y);
  520. }
  521. if( j < Gnd ) {
  522. gndgeom(j,rp,xp,yp);
  523. } else if( j >= Gnd && j < (Gnd + Rec) ) {
  524. recgeom(j,rp,xp,yp);
  525. } else if( j >= (Gnd + Rec) && j < (Gnd + Rec + Cir) ) {
  526. cirgeom(j,rp,xp,yp,alp);
  527. } else if( j >= (Gnd + Rec + Cir) && j < (Gnd + Rec + Cir + Tra) ) {
  528. trageom(j,rp,xp,yp);
  529. }
  530. giimn = gii(x,y,xp,yp);
  531. return(giimn);
  532. }
  533. // calculate the derivative of imaginary part of Gi
  534. double giipmn(int i, int j, double r, double rp) {
  535. double x,y,xp,yp,al,alp;
  536. double giipmn;
  537. int n,np;
  538. Vector tmpv(2);
  539. int ii = j-Gnd-Rec-Cir;
  540. if( i < Gnd ) {
  541. n = gndgeom(i,r,x,y);
  542. } else if( i >= Gnd && i < (Gnd + Rec) ) {
  543. n = recgeom(i,r,x,y);
  544. } else if( i >= (Gnd + Rec) && i < (Gnd + Rec + Cir) ) {
  545. n = cirgeom(i,r,x,y,al);
  546. } else if( i >= (Gnd + Rec + Cir) && i < (Gnd + Rec + Cir + Tra) ) {
  547. n = trageom(i,r,x,y);
  548. }
  549. if( j < Gnd ) {
  550. np = gndgeom(j,rp,xp,yp);
  551. } else if( j >= Gnd && j < (Gnd + Rec) ) {
  552. np = recgeom(j,rp,xp,yp);
  553. } else if( j >= (Gnd + Rec) && j < (Gnd + Rec + Cir) ) {
  554. np = cirgeom(j,rp,xp,yp,alp);
  555. } else if( j >= (Gnd + Rec + Cir) && j < (Gnd + Rec + Cir + Tra) ) {
  556. np = trageom(j,rp,xp,yp);
  557. }
  558. switch( np ) {
  559. case 1:
  560. giipmn = giiy(x,y,xp,yp);
  561. break;
  562. case 2:
  563. giipmn = -giix(x,y,xp,yp);
  564. break;
  565. case 3:
  566. giipmn = -giiy(x,y,xp,yp);
  567. break;
  568. case 4:
  569. giipmn = giix(x,y,xp,yp);
  570. break;
  571. case 5: // circle
  572. giipmn = -giix(x,y,xp,yp)*cos(alp) - giiy(x,y,xp,yp)*sin(alp);
  573. break;
  574. case 6:
  575. tmpv[0] = -giix(x,y,xp,yp);
  576. tmpv[1] = -giiy(x,y,xp,yp);
  577. giipmn = tmpv * tangle(Tradix(ii,0),Tradiy(ii,0),Tradix(ii,1),Tradiy(ii,1),Tradix(ii,2),Tradiy(ii,2));
  578. break;
  579. case 7:
  580. tmpv[0] = -giix(x,y,xp,yp);
  581. tmpv[1] = -giiy(x,y,xp,yp);
  582. giipmn = tmpv * tangle(Tradix(ii,1),Tradiy(ii,1),Tradix(ii,2),Tradiy(ii,2),Tradix(ii,3),Tradiy(ii,3));
  583. break;
  584. case 8:
  585. tmpv[0] = -giix(x,y,xp,yp);
  586. tmpv[1] = -giiy(x,y,xp,yp);
  587. giipmn = tmpv * tangle(Tradix(ii,2),Tradiy(ii,2),Tradix(ii,3),Tradiy(ii,3),Tradix(ii,0),Tradiy(ii,0));
  588. break;
  589. case 9: // 6-9 Trapezoid
  590. tmpv[0] = -giix(x,y,xp,yp);
  591. tmpv[1] = -giiy(x,y,xp,yp);
  592. giipmn = tmpv * tangle(Tradix(ii,3),Tradiy(ii,3),Tradix(ii,0),Tradiy(ii,0),Tradix(ii,1),Tradiy(ii,1));
  593. break;
  594. default:
  595. cerr << np << endl;
  596. cerr << "unidentified number" << endl;
  597. }
  598. return(giipmn);
  599. }
  600. // calculate the value of Gi,combine the real and imaginary parts
  601. Complex gimn(int i, int j, double r, double rp) {
  602. Complex gimn;
  603. gimn = girmn(i,j,r,rp)+Ic*giimn(i,j,r,rp);
  604. return(gimn);
  605. }
  606. // calculate the derivative of Gi, combine the real and imaginary parts.
  607. Complex gipmn(int i, int j, double r, double rp) {
  608. Complex gipmn;
  609. gipmn = girpmn(i,j,r,rp)+Ic*giipmn(i,j,r,rp);
  610. return(gipmn);
  611. }
  612. Vector tangle( double x1, double y1, double x2, double y2, double x3, double y3 ) {
  613. Vector v1(2), v2(2), v(2);
  614. v1[0] = x1-x2;
  615. v1[1] = y1-y2;
  616. v1 = v1/norm2(v1);
  617. v2[0] = x3-x2;
  618. v2[1] = y3-y2;
  619. v2 = v2/norm2(v2);
  620. v = (v1*v2)*v1 - v2;
  621. v = v/norm2(v);
  622. return(v);
  623. }