Startsida Utforska Hjälp
Logga in
dmalhotra
/
2024-06-13-talk-CSBQ
1
0
Fork 0
Filer Ärenden 0 Pull-förfrågningar 0 Wiki
Träd: b234a372d4
Grenar Taggar
2024-06-13-t...
/
code
/
BIEST
/
include
/
biest
/
surface.txx

surface.txx 23 KB
Historik Rå

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362 
  1. #include <biest/surface_op.hpp>
    2. 

  2. #include <sctl.hpp>
    3. 

  3. 

  4. namespace biest {
    5. 

  5. 

  6. template <class Real> Surface<Real>::Surface(sctl::Long Nt, sctl::Long Np, SurfType type) : Nt0(Nt), Np0(Np) {
    7. 

  7.   auto stell_geom = [](sctl::Vector<Real>& X0, sctl::Long Nt, sctl::Long Np, Real R0, Real a, Real b, bool ellipse) {
    8. 

  8.     X0.ReInit(COORD_DIM * Nt * Np);
    9. 

  9.     sctl::Matrix<Real> X(2,1), X_(2,1);
    10. 

  10.     sctl::Matrix<Real> R(2,2), Rt(2,2);
    11. 

  11.     for (sctl::Long t = 0; t < Nt; t++) {
    12. 

  12.       R[0][0] = sctl::cos(1.5*2*sctl::const_pi<Real>()*t/Nt);
    13. 

  13.       R[0][1] = sctl::sin(1.5*2*sctl::const_pi<Real>()*t/Nt);
    14. 

  14.       R[1][0] =-sctl::sin(1.5*2*sctl::const_pi<Real>()*t/Nt);
    15. 

  15.       R[1][1] = sctl::cos(1.5*2*sctl::const_pi<Real>()*t/Nt);
    16. 

  16. 

  17.       Rt[0][0] = sctl::cos(1.5*2*sctl::const_pi<Real>()*t/Nt);
    18. 

  18.       Rt[1][0] = sctl::sin(1.5*2*sctl::const_pi<Real>()*t/Nt);
    19. 

  19.       Rt[0][1] =-sctl::sin(1.5*2*sctl::const_pi<Real>()*t/Nt);
    20. 

  20.       Rt[1][1] = sctl::cos(1.5*2*sctl::const_pi<Real>()*t/Nt);
    21. 

  21.       for (sctl::Long p = 0; p < Np; p++) {
    22. 

  22.         X[0][0] = sctl::cos(2 * sctl::const_pi<Real>() * p / Np);
    23. 

  23.         X[1][0] = sctl::sin(2 * sctl::const_pi<Real>() * p / Np);
    24. 

  24.         if (!ellipse) X[0][0] = (X[0][0] + sctl::exp(sctl::cos(2 * sctl::const_pi<Real>() * p / Np)) - 1.5) / 2.4;
    25. 

  25. 

  26.         sctl::Matrix<Real>::GEMM(X_, R, X);
    27. 

  27.         X_[0][0] *= a;
    28. 

  28.         X_[1][0] *= b;
    29. 

  29.         sctl::Matrix<Real>::GEMM(X, Rt, X_);
    30. 

  30. 

  31.         X0[(0 * Nt + t) * Np + p] = (X[0][0] + (Real)R0) * sctl::cos(2 * sctl::const_pi<Real>() * t / Nt);
    32. 

  32.         X0[(1 * Nt + t) * Np + p] = (X[0][0] + (Real)R0) * sctl::sin(2 * sctl::const_pi<Real>() * t / Nt);
    33. 

  33.         X0[(2 * Nt + t) * Np + p] = (X[1][0] + (Real) 0);
    34. 

  34.       }
    35. 

  35.     }
    36. 

  36.   };
    37. 

  37.   if (type == SurfType::AxisymCircleWide) {
    38. 

  38.     stell_geom(X0_, Nt0, Np0, 2.0, 1.0, 1.0, true);
    39. 

  39.   } else if (type == SurfType::AxisymCircleNarrow) {
    40. 

  40.     stell_geom(X0_, Nt0, Np0, 2.0, 0.5, 0.5, true);
    41. 

  41.   } else if (type == SurfType::AxisymWide) {
    42. 

  42.     stell_geom(X0_, Nt0, Np0, 2.0, 1.0, 1.0, false);
    43. 

  43.   } else if (type == SurfType::AxisymNarrow) {
    44. 

  44.     stell_geom(X0_, Nt0, Np0, 2.0, 0.5, 0.5, false);
    45. 

  45.   } else if (type == SurfType::RotatingEllipseWide) {
    46. 

  46.     stell_geom(X0_, Nt0, Np0, 2.0, 0.7, 1.0, true);
    47. 

  47.   } else if (type == SurfType::RotatingEllipseNarrow) {
    48. 

  48.     stell_geom(X0_, Nt0, Np0, 2.0, 0.3, 0.55, true);
    49. 

  49.   } else if (type == SurfType::Quas3) {
    50. 

  50.     sctl::Matrix<double> X, Y, Z;
    51. 

  51.     X.Read("geom/Quas3-X.mat");
    52. 

  52.     Y.Read("geom/Quas3-Y.mat");
    53. 

  53.     Z.Read("geom/Quas3-Z.mat");
    54. 

  54.     sctl::Long Nt = X.Dim(0)-1, Np = X.Dim(1)-1;
    55. 

  55.     sctl::Vector<Real> XX(COORD_DIM * Nt * Np);
    56. 

  56.     for (sctl::Long t = 0; t < Nt; t++) {
    57. 

  57.       for (sctl::Long p = 0; p < Np; p++) {
    58. 

  58.         XX[(0 * Nt + t) * Np + p] = (Real)X[t][p] * 0.45;
    59. 

  59.         XX[(1 * Nt + t) * Np + p] = (Real)Y[t][p] * 0.45;
    60. 

  60.         XX[(2 * Nt + t) * Np + p] = (Real)Z[t][p] * 0.45;
    61. 

  61.       }
    62. 

  62.     }
    63. 

  63.     SurfaceOp<Real>::Upsample(XX, Nt, Np, this->Coord(), Nt0, Np0);
    64. 

  64.   } else if (type == SurfType::LHD) {
    65. 

  65.     sctl::Matrix<double> X, Y, Z;
    66. 

  66.     X.Read("geom/LHD-X.mat");
    67. 

  67.     Y.Read("geom/LHD-Y.mat");
    68. 

  68.     Z.Read("geom/LHD-Z.mat");
    69. 

  69.     sctl::Long Nt = X.Dim(0), Np = X.Dim(1);
    70. 

  70.     sctl::Vector<Real> XX(COORD_DIM * Nt * Np);
    71. 

  71.     for (sctl::Long t = 0; t < Nt; t++) {
    72. 

  72.       for (sctl::Long p = 0; p < Np; p++) {
    73. 

  73.         XX[(0 * Nt + t) * Np + p] = (Real)X[t][p] * 0.25;
    74. 

  74.         XX[(1 * Nt + t) * Np + p] = (Real)Y[t][p] * 0.25;
    75. 

  75.         XX[(2 * Nt + t) * Np + p] = (Real)Z[t][p] * 0.25;
    76. 

  76.       }
    77. 

  77.     }
    78. 

  78.     SurfaceOp<Real>::Upsample(XX, Nt, Np, this->Coord(), Nt0, Np0);
    79. 

  79.   } else if (type == SurfType::W7X) {
    80. 

  80.     sctl::Matrix<double> X, Y, Z;
    81. 

  81.     X.Read("geom/W7X-X.mat");
    82. 

  82.     Y.Read("geom/W7X-Y.mat");
    83. 

  83.     Z.Read("geom/W7X-Z.mat");
    84. 

  84.     sctl::Long Nt = X.Dim(0), Np = X.Dim(1);
    85. 

  85.     sctl::Vector<Real> XX(COORD_DIM * Nt * Np);
    86. 

  86.     for (sctl::Long t = 0; t < Nt; t++) {
    87. 

  87.       for (sctl::Long p = 0; p < Np; p++) {
    88. 

  88.         XX[(0 * Nt + t) * Np + p] = (Real)X[t][p] * 0.45;
    89. 

  89.         XX[(1 * Nt + t) * Np + p] = (Real)Y[t][p] * 0.45;
    90. 

  90.         XX[(2 * Nt + t) * Np + p] = (Real)Z[t][p] * 0.45;
    91. 

  91.       }
    92. 

  92.     }
    93. 

  93.     SurfaceOp<Real>::Upsample(XX, Nt, Np, this->Coord(), Nt0, Np0);
    94. 

  94.   } else if (type == SurfType::Stell) {
    95. 

  95.     sctl::Matrix<Real> Rcoeff(2,2), Zcoeff(2,2);
    96. 

  96.     Rcoeff = 0;
    97. 

  97.     Zcoeff = 0;
    98. 

  98.     Rcoeff[0][0] =  1.0e+01;
    99. 

  99.     Rcoeff[0][1] =  1.0e+00;
    100. 

  100.     Rcoeff[1][1] =  2.5e-01;
    101. 

  101.     Zcoeff[0][1] = -1.0e+00;
    102. 

  102.     Zcoeff[1][1] =  2.5e-01;
    103. 

  103. 

  104.     X0_.ReInit(COORD_DIM * Nt0 * Np0);
    105. 

  105.     for (long t = 0; t < Nt0; t++) {
    106. 

  106.       for (long p = 0; p < Np0; p++) {
    107. 

  107.         Real theta = t * 2 * sctl::const_pi<Real>() / Nt0;
    108. 

  108.         Real phi   = p * 2 * sctl::const_pi<Real>() / Np0;
    109. 

  109.         Real R = 0, Z = 0;
    110. 

  110.         for (long i = 0; i < 2; i++) {
    111. 

  111.           for (long j = 0; j < 2; j++) {
    112. 

  112.             R += Rcoeff[i][j] * sctl::cos(j*phi - 5*i*theta);
    113. 

  113.             Z += Zcoeff[i][j] * sctl::sin(j*phi - 5*i*theta);
    114. 

  114.           }
    115. 

  115.         }
    116. 

  116.         X0_[0*Nt0*Np0 + t*Np0 + p] = R * sctl::cos(theta);
    117. 

  117.         X0_[1*Nt0*Np0 + t*Np0 + p] = R * sctl::sin(theta);
    118. 

  118.         X0_[2*Nt0*Np0 + t*Np0 + p] = Z;
    119. 

  119.       }
    120. 

  120.     }
    121. 

  121.   } else if (type == SurfType::W7X_) {
    122. 

  122.     sctl::Matrix<Real> Rcoeff(21,21), Zcoeff(21,21);
    123. 

  123.     Rcoeff = 0;
    124. 

  124.     Zcoeff = 0;
    125. 

  125.     auto Rbc = [&Rcoeff](sctl::Long i, sctl::Long j) ->Real& { return Rcoeff[i+10][j+10]; };
    126. 

  126.     auto Zbs = [&Zcoeff](sctl::Long i, sctl::Long j) ->Real& { return Zcoeff[i+10][j+10]; };
    127. 

  127.     { // Set Rcoeff, Zcoeff
    128. 

  128.       Rbc(0,0)    =  5.499222200000000E+00; Zbs(0,0)    =  0.000000000000000E+00;
    129. 

  129.       Rbc(1,0)    =  2.862415760000000E-01; Zbs(1,0)    = -2.476328720000000E-01;
    130. 

  130.       Rbc(2,0)    = -2.821010320000000E-03; Zbs(2,0)    = -3.798656180000000E-03;
    131. 

  131.       Rbc(3,0)    =  4.622923140000000E-04; Zbs(3,0)    =  1.379383388000000E-03;
    132. 

  132.       Rbc(4,0)    = -6.360328500000000E-04; Zbs(4,0)    =  6.992897620000000E-04;
    133. 

  133.       Rbc(5,0)    = -1.135627524000000E-04; Zbs(5,0)    = -3.409457720000000E-05;
    134. 

  134.       Rbc(6,0)    = -1.816774936000000E-04; Zbs(6,0)    = -1.255490184000000E-04;
    135. 

  135.       Rbc(7,0)    =  1.341661870000000E-04; Zbs(7,0)    =  1.037604320000000E-04;
    136. 

  136.       Rbc(8,0)    = -2.663143840000000E-04; Zbs(8,0)    = -4.829788340000000E-04;
    137. 

  137.       Rbc(9,0)    =  7.395449780000000E-05; Zbs(9,0)    =  2.269942520000000E-04;
    138. 

  138.       Rbc(10,0)   =  1.213701750000000E-04; Zbs(10,0)   =  1.063213162000000E-04;
    139. 

  139.       Rbc(-10,1)  =  3.149812680000000E-05; Zbs(-10,1)  =  5.847707230000000E-05;
    140. 

  140.       Rbc(-8,1)   = -3.424627190000000E-04; Zbs(-8,1)   =  2.389303750000000E-04;
    141. 

  141.       Rbc(-7,1)   = -5.528744530000000E-05; Zbs(-7,1)   = -1.910728510000000E-04;
    142. 

  142.       Rbc(-6,1)   =  6.117681689999999E-05; Zbs(-6,1)   = -5.092602740000000E-04;
    143. 

  143.       Rbc(-5,1)   =  3.696634690000000E-04; Zbs(-5,1)   = -3.176068510000000E-04;
    144. 

  144.       Rbc(-4,1)   =  2.191307230000000E-04; Zbs(-4,1)   =  7.293541070000000E-04;
    145. 

  145.       Rbc(-3,1)   =  1.269562650000000E-03; Zbs(-3,1)   = -6.193015160000001E-04;
    146. 

  146.       Rbc(-2,1)   = -2.148173690000000E-02; Zbs(-2,1)   = -2.227519550000000E-02;
    147. 

  147.       Rbc(-1,1)   = -2.080619220000000E-01; Zbs(-1,1)   = -1.793381930000000E-01;
    148. 

  148.       Rbc(0,1)    =  4.591451160000000E-01; Zbs(0,1)    = -6.108126480000000E-01;
    149. 

  149.       Rbc(1,1)    =  3.971068830000000E-02; Zbs(1,1)    = -3.707128920000000E-02;
    150. 

  150.       Rbc(2,1)    =  3.618299070000000E-03; Zbs(2,1)    = -5.337647490000000E-03;
    151. 

  151.       Rbc(3,1)    = -3.891312730000000E-04; Zbs(3,1)    = -3.850478840000000E-04;
    152. 

  152.       Rbc(4,1)    =  4.101864840000000E-04; Zbs(4,1)    = -5.420620840000000E-04;
    153. 

  153.       Rbc(5,1)    =  3.192006150000000E-05; Zbs(5,1)    =  5.585418590000000E-04;
    154. 

  154.       Rbc(6,1)    = -5.386882030000000E-04; Zbs(6,1)    = -2.760504880000000E-04;
    155. 

  155.       Rbc(7,1)    = -3.173680500000000E-04; Zbs(7,1)    =  7.432511449999999E-05;
    156. 

  156.       Rbc(8,1)    =  4.024401190000000E-04; Zbs(8,1)    = -1.269831920000000E-04;
    157. 

  157.       Rbc(9,1)    =  2.130320640000000E-04; Zbs(9,1)    = -1.203549150000000E-04;
    158. 

  158.       Rbc(10,1)   =  1.636947630000000E-04; Zbs(10,1)   =  8.987937330000000E-05;
    159. 

  159.       Rbc(-10,2)  = -7.765489310000000E-05; Zbs(-10,2)  = -8.678722490000000E-05;
    160. 

  160.       Rbc(-8,2)   =  4.456479760000000E-05; Zbs(-8,2)   =  1.145912760000000E-04;
    161. 

  161.       Rbc(-7,2)   =  5.600939840000000E-05; Zbs(-7,2)   = -2.025932200000000E-06;
    162. 

  162.       Rbc(-6,2)   =  3.165439980000000E-05; Zbs(-6,2)   =  3.892001780000000E-05;
    163. 

  163.       Rbc(-5,2)   =  2.876781000000000E-04; Zbs(-5,2)   =  3.010497420000000E-04;
    164. 

  164.       Rbc(-4,2)   = -3.334248370000000E-04; Zbs(-4,2)   =  2.345406280000000E-04;
    165. 

  165.       Rbc(-3,2)   =  5.907835000000000E-04; Zbs(-3,2)   =  2.120520070000000E-04;
    166. 

  166.       Rbc(-2,2)   =  5.679743510000000E-02; Zbs(-2,2)   =  5.802810780000000E-02;
    167. 

  167.       Rbc(-1,2)   =  3.730642040000000E-02; Zbs(-1,2)   = -2.576171930000000E-02;
    168. 

  168.       Rbc(0,2)    =  5.279093270000000E-02; Zbs(0,2)    = -1.753837990000000E-02;
    169. 

  169.       Rbc(1,2)    =  1.057394420000000E-02; Zbs(1,2)    = -5.454225200000000E-03;
    170. 

  170.       Rbc(2,2)    =  2.688744090000000E-03; Zbs(2,2)    = -2.870790520000000E-04;
    171. 

  171.       Rbc(3,2)    =  8.698832710000000E-04; Zbs(3,2)    = -6.835677030000000E-04;
    172. 

  172.       Rbc(4,2)    =  1.170861370000000E-04; Zbs(4,2)    = -2.436806940000000E-04;
    173. 

  173.       Rbc(5,2)    = -3.468535580000000E-04; Zbs(5,2)    = -3.554177930000000E-04;
    174. 

  174.       Rbc(6,2)    =  3.854471210000000E-04; Zbs(6,2)    =  4.853309450000000E-04;
    175. 

  175.       Rbc(7,2)    = -1.035513210000000E-04; Zbs(7,2)    = -3.634096180000000E-04;
    176. 

  176.       Rbc(8,2)    = -1.028097890000000E-04; Zbs(8,2)    =  1.530850250000000E-04;
    177. 

  177.       Rbc(9,2)    =  1.632011900000000E-04; Zbs(9,2)    = -1.707748690000000E-04;
    178. 

  178.       Rbc(10,2)   =  1.581787490000000E-04; Zbs(10,2)   = -3.585159010000000E-05;
    179. 

  179.       Rbc(-10,3)  =  4.771322000000000E-05; Zbs(-10,3)  =  9.188864070000000E-05;
    180. 

  180.       Rbc(-8,3)   = -6.313692220000000E-05; Zbs(-8,3)   = -2.572672720000000E-05;
    181. 

  181.       Rbc(-7,3)   =  5.097030830000000E-05; Zbs(-7,3)   =  3.932101880000000E-05;
    182. 

  182.       Rbc(-6,3)   = -9.952454240000000E-05; Zbs(-6,3)   = -9.632583980000000E-05;
    183. 

  183.       Rbc(-5,3)   = -2.012500210000000E-04; Zbs(-5,3)   = -1.439074920000000E-04;
    184. 

  184.       Rbc(-4,3)   = -5.638133750000001E-04; Zbs(-4,3)   = -3.602899270000000E-04;
    185. 

  185.       Rbc(-3,3)   = -1.249071190000000E-02; Zbs(-3,3)   = -1.229819730000000E-02;
    186. 

  186.       Rbc(-2,3)   = -1.470456180000000E-02; Zbs(-2,3)   = -1.098330420000000E-02;
    187. 

  187.       Rbc(-1,3)   = -2.963558930000000E-03; Zbs(-1,3)   = -3.655102520000000E-03;
    188. 

  188.       Rbc(0,3)    =  3.531728600000000E-03; Zbs(0,3)    = -1.279854780000000E-03;
    189. 

  189.       Rbc(1,3)    =  1.401374340000000E-03; Zbs(1,3)    =  1.590869290000000E-03;
    190. 

  190.       Rbc(2,3)    = -1.273655750000000E-04; Zbs(2,3)    = -1.929131170000000E-04;
    191. 

  191.       Rbc(3,3)    =  7.079597950000001E-05; Zbs(3,3)    = -1.224304200000000E-04;
    192. 

  192.       Rbc(4,3)    =  2.850121260000000E-04; Zbs(4,3)    =  2.598242880000000E-04;
    193. 

  193.       Rbc(5,3)    =  1.999424460000000E-04; Zbs(5,3)    =  1.123146280000000E-04;
    194. 

  194.       Rbc(6,3)    = -3.092845690000000E-04; Zbs(6,3)    = -2.016395580000000E-04;
    195. 

  195.       Rbc(7,3)    = -5.220858140000000E-05; Zbs(7,3)    =  1.316453710000000E-04;
    196. 

  196.       Rbc(8,3)    =  2.440456610000000E-04; Zbs(8,3)    =  3.632665850000000E-05;
    197. 

  197.       Rbc(9,3)    =  1.636232450000000E-04; Zbs(9,3)    = -6.465022640000000E-05;
    198. 

  198.       Rbc(10,3)   =  6.580775879999999E-05; Zbs(10,3)   = -5.470819060000000E-05;
    199. 

  199.       Rbc(-10,4)  = -1.584638410000000E-06; Zbs(-10,4)  =  6.451981970000000E-06;
    200. 

  200.       Rbc(-8,4)   = -1.656855160000000E-05; Zbs(-8,4)   = -7.496657650000000E-05;
    201. 

  201.       Rbc(-7,4)   =  1.657078710000000E-05; Zbs(-7,4)   =  2.015102120000000E-05;
    202. 

  202.       Rbc(-6,4)   =  1.152800160000000E-04; Zbs(-6,4)   =  5.273493610000000E-05;
    203. 

  203.       Rbc(-5,4)   =  1.587317930000000E-04; Zbs(-5,4)   = -1.770118060000000E-06;
    204. 

  204.       Rbc(-4,4)   =  1.693684220000000E-03; Zbs(-4,4)   =  1.459908040000000E-03;
    205. 

  205.       Rbc(-3,4)   =  5.456882450000000E-03; Zbs(-3,4)   =  5.427735220000000E-03;
    206. 

  206.       Rbc(-2,4)   =  1.225123000000000E-03; Zbs(-2,4)   = -5.642061490000000E-04;
    207. 

  207.       Rbc(-1,4)   = -1.634672880000000E-03; Zbs(-1,4)   = -1.759500600000000E-03;
    208. 

  208.       Rbc(0,4)    =  1.972846670000000E-03; Zbs(0,4)    = -9.862657130000000E-04;
    209. 

  209.       Rbc(1,4)    = -3.061925850000000E-04; Zbs(1,4)    =  6.134663420000001E-04;
    210. 

  210.       Rbc(2,4)    = -2.248673790000000E-04; Zbs(2,4)    =  9.844427600000000E-05;
    211. 

  211.       Rbc(3,4)    =  1.440742170000000E-04; Zbs(3,4)    =  2.386755350000000E-04;
    212. 

  212.       Rbc(4,4)    = -3.970536280000000E-05; Zbs(4,4)    = -5.037994770000000E-05;
    213. 

  213.       Rbc(5,4)    = -1.067147990000000E-04; Zbs(5,4)    =  5.683708790000000E-05;
    214. 

  214.       Rbc(6,4)    =  8.551267700000000E-05; Zbs(6,4)    =  5.857943570000000E-05;
    215. 

  215.       Rbc(7,4)    = -1.050568620000000E-04; Zbs(7,4)    = -1.372536220000000E-04;
    216. 

  216.       Rbc(8,4)    = -5.167821720000000E-05; Zbs(8,4)    = -7.182784940000000E-05;
    217. 

  217.       Rbc(9,4)    =  2.259193970000000E-06; Zbs(9,4)    = -2.070346460000000E-05;
    218. 

  218.       Rbc(10,4)   =  2.383410400000000E-06; Zbs(10,4)   = -1.711598370000000E-05;
    219. 

  219.       Rbc(-10,5)  = -1.054277140000000E-05; Zbs(-10,5)  = -4.590453350000000E-05;
    220. 

  220.       Rbc(-8,5)   =  1.840472560000000E-05; Zbs(-8,5)   =  2.408429910000000E-05;
    221. 

  221.       Rbc(-7,5)   = -2.421259200000000E-05; Zbs(-7,5)   = -3.838815110000000E-05;
    222. 

  222.       Rbc(-6,5)   = -7.576990619999999E-05; Zbs(-6,5)   = -3.713491020000000E-05;
    223. 

  223.       Rbc(-5,5)   = -2.754279550000000E-04; Zbs(-5,5)   = -2.168078010000000E-04;
    224. 

  224.       Rbc(-4,5)   = -1.706462950000000E-03; Zbs(-4,5)   = -1.644204180000000E-03;
    225. 

  225.       Rbc(-3,5)   = -1.313493520000000E-03; Zbs(-3,5)   = -9.164101610000000E-04;
    226. 

  226.       Rbc(-2,5)   =  1.549990410000000E-04; Zbs(-2,5)   =  3.455799460000000E-04;
    227. 

  227.       Rbc(-1,5)   = -1.044131360000000E-03; Zbs(-1,5)   = -4.688486220000000E-05;
    228. 

  228.       Rbc(0,5)    =  6.811915640000000E-04; Zbs(0,5)    = -5.596368950000000E-04;
    229. 

  229.       Rbc(1,5)    = -2.239129020000000E-04; Zbs(1,5)    =  1.351560670000000E-04;
    230. 

  230.       Rbc(2,5)    = -9.469807639999999E-05; Zbs(2,5)    =  1.702559280000000E-05;
    231. 

  231.       Rbc(3,5)    = -2.114701210000000E-05; Zbs(3,5)    = -3.835918030000000E-06;
    232. 

  232.       Rbc(4,5)    =  1.354520640000000E-04; Zbs(4,5)    =  1.413910340000000E-04;
    233. 

  233.       Rbc(5,5)    = -1.082798510000000E-04; Zbs(5,5)    = -1.868524930000000E-04;
    234. 

  234.       Rbc(6,5)    = -4.385285270000000E-05; Zbs(6,5)    =  3.311395990000000E-06;
    235. 

  235.       Rbc(7,5)    = -1.993355920000000E-05; Zbs(7,5)    = -2.425407170000000E-05;
    236. 

  236.       Rbc(8,5)    = -1.381184910000000E-05; Zbs(8,5)    = -3.354113410000000E-05;
    237. 

  237.       Rbc(9,5)    = -8.047951870000001E-06; Zbs(9,5)    = -2.617354950000000E-05;
    238. 

  238.       Rbc(10,5)   =  5.442943530000000E-06; Zbs(10,5)   = -7.540077230000000E-06;
    239. 

  239.       Rbc(-10,6)  =  5.897904990000000E-06; Zbs(-10,6)  =  1.015941240000000E-05;
    240. 

  240.       Rbc(-8,6)   = -1.711315030000000E-06; Zbs(-8,6)   =  1.460063060000000E-05;
    241. 

  241.       Rbc(-7,6)   =  2.779279110000000E-05; Zbs(-7,6)   =  9.346117740000001E-06;
    242. 

  242.       Rbc(-6,6)   =  6.159498210000000E-05; Zbs(-6,6)   =  3.757257530000000E-05;
    243. 

  243.       Rbc(-5,6)   =  3.838807420000000E-04; Zbs(-5,6)   =  3.538485060000000E-04;
    244. 

  244.       Rbc(-4,6)   =  5.812892710000000E-04; Zbs(-4,6)   =  5.854490200000000E-04;
    245. 

  245.       Rbc(-3,6)   = -1.282036420000000E-05; Zbs(-3,6)   = -9.995316180000000E-05;
    246. 

  246.       Rbc(-2,6)   =  3.747466590000000E-04; Zbs(-2,6)   =  4.378838950000000E-04;
    247. 

  247.       Rbc(-1,6)   =  4.290153220000000E-04; Zbs(-1,6)   = -7.532518330000000E-04;
    248. 

  248.       Rbc(0,6)    =  8.328937110000000E-04; Zbs(0,6)    = -6.075503770000000E-04;
    249. 

  249.       Rbc(1,6)    = -1.772339600000000E-05; Zbs(1,6)    = -7.259493060000000E-05;
    250. 

  250.       Rbc(2,6)    =  3.629202880000000E-05; Zbs(2,6)    =  5.097636310000000E-05;
    251. 

  251.       Rbc(3,6)    =  9.744514720000001E-06; Zbs(3,6)    = -2.890839380000000E-05;
    252. 

  252.       Rbc(4,6)    =  3.762214890000000E-06; Zbs(4,6)    = -2.016824820000000E-05;
    253. 

  253.       Rbc(5,6)    = -2.473214420000000E-05; Zbs(5,6)    =  4.460497530000000E-06;
    254. 

  254.       Rbc(6,6)    =  9.393233570000001E-06; Zbs(6,6)    = -1.180516910000000E-05;
    255. 

  255.       Rbc(7,6)    = -7.516850280000000E-06; Zbs(7,6)    = -1.452989350000000E-05;
    256. 

  256.       Rbc(8,6)    = -1.256492920000000E-06; Zbs(8,6)    = -1.179898990000000E-05;
    257. 

  257.       Rbc(9,6)    = -3.141617150000000E-06; Zbs(9,6)    = -1.756033910000000E-05;
    258. 

  258.       Rbc(10,6)   = -2.246448250000000E-06; Zbs(10,6)   = -9.158244980000000E-06;
    259. 

  259.       Rbc(-10,7)  =  5.929060270000000E-06; Zbs(-10,7)  = -1.076521400000000E-05;
    260. 

  260.       Rbc(-8,7)   = -4.097376180000000E-06; Zbs(-8,7)   = -7.849626889999999E-06;
    261. 

  261.       Rbc(-7,7)   = -1.789842870000000E-05; Zbs(-7,7)   = -5.481835500000000E-06;
    262. 

  262.       Rbc(-6,7)   = -7.895979130000000E-05; Zbs(-6,7)   = -6.009604550000000E-05;
    263. 

  263.       Rbc(-5,7)   = -1.699260490000000E-04; Zbs(-5,7)   = -2.072394130000000E-04;
    264. 

  264.       Rbc(-4,7)   = -6.638575870000000E-05; Zbs(-4,7)   = -3.656338330000000E-05;
    265. 

  265.       Rbc(-3,7)   =  9.249536449999999E-05; Zbs(-3,7)   =  1.022424950000000E-04;
    266. 

  266.       Rbc(-2,7)   =  1.172203510000000E-03; Zbs(-2,7)   =  1.077845730000000E-03;
    267. 

  267.       Rbc(-1,7)   = -2.573141910000000E-03; Zbs(-1,7)   =  3.037810800000000E-03;
    268. 

  268.       Rbc(0,7)    =  6.449810380000000E-05; Zbs(0,7)    = -3.584656230000000E-05;
    269. 

  269.       Rbc(1,7)    =  2.190318120000000E-05; Zbs(1,7)    = -2.286948060000000E-05;
    270. 

  270.       Rbc(2,7)    =  1.722638980000000E-05; Zbs(2,7)    = -1.325952420000000E-05;
    271. 

  271.       Rbc(3,7)    = -2.541553570000000E-05; Zbs(3,7)    =  4.086201040000000E-06;
    272. 

  272.       Rbc(4,7)    =  1.336724870000000E-05; Zbs(4,7)    = -2.288792230000000E-05;
    273. 

  273.       Rbc(5,7)    = -7.251710020000000E-06; Zbs(5,7)    = -1.480782940000000E-05;
    274. 

  274.       Rbc(6,7)    =  7.086029940000000E-06; Zbs(6,7)    =  3.959755080000000E-07;
    275. 

  275.       Rbc(7,7)    = -2.160257930000000E-05; Zbs(7,7)    =  1.387646540000000E-05;
    276. 

  276.       Rbc(8,7)    = -7.244894970000000E-06; Zbs(8,7)    =  5.399568120000000E-06;
    277. 

  277.       Rbc(9,7)    = -3.206902820000000E-06; Zbs(9,7)    = -6.539029700000000E-06;
    278. 

  278.       Rbc(10,7)   =  1.060093110000000E-05; Zbs(10,7)   = -9.468079560000000E-06;
    279. 

  279.       Rbc(-10,8)  = -1.145540760000000E-06; Zbs(-10,8)  = -1.426725280000000E-05;
    280. 

  280.       Rbc(-8,8)   =  8.075915779999999E-07; Zbs(-8,8)   = -2.090075330000000E-06;
    281. 

  281.       Rbc(-7,8)   =  1.225482070000000E-05; Zbs(-7,8)   =  1.002867410000000E-05;
    282. 

  282.       Rbc(-6,8)   =  4.185049880000000E-05; Zbs(-6,8)   =  3.627089890000000E-05;
    283. 

  283.       Rbc(-5,8)   =  5.202132510000000E-05; Zbs(-5,8)   =  3.298637690000000E-05;
    284. 

  284.       Rbc(-4,8)   = -1.790487510000000E-05; Zbs(-4,8)   = -2.696253960000000E-05;
    285. 

  285.       Rbc(-3,8)   = -6.164140990000000E-04; Zbs(-3,8)   = -6.284073990000000E-04;
    286. 

  286.       Rbc(-2,8)   = -3.576911990000000E-04; Zbs(-2,8)   =  3.465212410000000E-04;
    287. 

  287.       Rbc(-1,8)   = -5.728932250000000E-04; Zbs(-1,8)   =  3.377205970000000E-04;
    288. 

  288.       Rbc(0,8)    = -8.260321650000000E-05; Zbs(0,8)    =  5.832294860000000E-05;
    289. 

  289.       Rbc(1,8)    = -3.522536740000000E-05; Zbs(1,8)    =  2.406038780000000E-05;
    290. 

  290.       Rbc(2,8)    =  2.684092340000000E-07; Zbs(2,8)    =  1.991160860000000E-05;
    291. 

  291.       Rbc(3,8)    =  1.192629450000000E-06; Zbs(3,8)    =  1.025583730000000E-07;
    292. 

  292.       Rbc(4,8)    =  1.060654930000000E-05; Zbs(4,8)    =  8.367961500000000E-06;
    293. 

  293.       Rbc(5,8)    = -1.069578210000000E-05; Zbs(5,8)    = -1.167514120000000E-05;
    294. 

  294.       Rbc(6,8)    =  5.935593950000000E-07; Zbs(6,8)    =  6.173553520000000E-06;
    295. 

  295.       Rbc(7,8)    = -1.782560770000000E-06; Zbs(7,8)    =  4.560812670000000E-07;
    296. 

  296.       Rbc(8,8)    = -2.019758340000000E-06; Zbs(8,8)    =  5.240396660000000E-06;
    297. 

  297.       Rbc(9,8)    = -2.664745030000000E-06; Zbs(9,8)    = -2.290649840000000E-06;
    298. 

  298.       Rbc(10,8)   =  2.187762770000000E-06; Zbs(10,8)   =  2.349449310000000E-06;
    299. 

  299.       Rbc(-10,9)  =  1.072796270000000E-06; Zbs(-10,9)  =  7.710368280000000E-07;
    300. 

  300.       Rbc(-8,9)   = -3.954187380000000E-06; Zbs(-8,9)   = -1.287612760000000E-06;
    301. 

  301.       Rbc(-7,9)   = -5.975561270000000E-06; Zbs(-7,9)   = -3.802088860000000E-06;
    302. 

  302.       Rbc(-6,9)   = -9.711522580000000E-06; Zbs(-6,9)   = -1.111428870000000E-05;
    303. 

  303.       Rbc(-5,9)   =  1.283299810000000E-05; Zbs(-5,9)   =  1.270473750000000E-05;
    304. 

  304.       Rbc(-4,9)   =  2.050599600000000E-04; Zbs(-4,9)   =  2.090269370000000E-04;
    305. 

  305.       Rbc(-3,9)   =  2.150777160000000E-04; Zbs(-3,9)   =  1.663527070000000E-04;
    306. 

  306.       Rbc(-2,9)   =  4.308369360000000E-05; Zbs(-2,9)   =  1.255918080000000E-04;
    307. 

  307.       Rbc(-1,9)   = -4.735076660000000E-05; Zbs(-1,9)   =  2.779910190000000E-06;
    308. 

  308.       Rbc(0,9)    =  1.787987800000000E-05; Zbs(0,9)    = -1.003487210000000E-04;
    309. 

  309.       Rbc(1,9)    = -3.850292090000000E-06; Zbs(1,9)    = -1.737272460000000E-06;
    310. 

  310.       Rbc(2,9)    = -6.038516310000000E-06; Zbs(2,9)    = -1.027369980000000E-05;
    311. 

  311.       Rbc(3,9)    = -5.018549580000000E-06; Zbs(3,9)    = -1.159928700000000E-05;
    312. 

  312.       Rbc(4,9)    =  4.292816850000000E-06; Zbs(4,9)    = -3.788971130000000E-06;
    313. 

  313.       Rbc(5,9)    =  6.944764760000000E-06; Zbs(5,9)    = -5.514237360000000E-07;
    314. 

  314.       Rbc(6,9)    =  4.340471530000000E-06; Zbs(6,9)    = -5.532296800000000E-06;
    315. 

  315.       Rbc(7,9)    = -5.070625820000000E-06; Zbs(7,9)    = -4.607223520000000E-06;
    316. 

  316.       Rbc(8,9)    = -2.124002090000000E-06; Zbs(8,9)    = -3.273312010000000E-06;
    317. 

  317.       Rbc(9,9)    = -2.177449410000000E-07; Zbs(9,9)    = -5.030283200000000E-06;
    318. 

  318.       Rbc(10,9)   = -4.607967820000000E-07; Zbs(10,9)   = -3.725635690000000E-06;
    319. 

  319.       Rbc(-10,10) = -1.281240510000000E-06; Zbs(-10,10) = -6.404339920000000E-06;
    320. 

  320.       Rbc(-8,10)  =  1.956168450000000E-06; Zbs(-8,10)  = -1.805927600000000E-06;
    321. 

  321.       Rbc(-7,10)  =  1.101080660000000E-06; Zbs(-7,10)  = -3.471541870000000E-07;
    322. 

  322.       Rbc(-6,10)  = -3.993429860000000E-06; Zbs(-6,10)  = -4.649932760000000E-06;
    323. 

  323.       Rbc(-5,10)  = -3.274758810000000E-05; Zbs(-5,10)  = -3.510246990000000E-05;
    324. 

  324.       Rbc(-4,10)  = -1.059813770000000E-04; Zbs(-4,10)  = -1.105165310000000E-04;
    325. 

  325.       Rbc(-3,10)  =  1.095947770000000E-05; Zbs(-3,10)  =  4.305883900000000E-05;
    326. 

  326.       Rbc(-2,10)  =  2.631080330000000E-04; Zbs(-2,10)  =  2.624212580000000E-04;
    327. 

  327.       Rbc(-1,10)  = -1.383145710000000E-04; Zbs(-1,10)  =  1.209196340000000E-04;
    328. 

  328.       Rbc(0,10)   =  1.047045790000000E-05; Zbs(0,10)   = -2.287772640000000E-05;
    329. 

  329.       Rbc(1,10)   =  2.905086960000000E-06; Zbs(1,10)   =  1.019706800000000E-06;
    330. 

  330.       Rbc(2,10)   = -7.151778040000000E-07; Zbs(2,10)   = -3.808867920000000E-06;
    331. 

  331.       Rbc(3,10)   = -2.321896650000000E-06; Zbs(3,10)   = -7.137445440000000E-07;
    332. 

  332.       Rbc(4,10)   =  1.611615800000000E-06; Zbs(4,10)   = -6.657593750000000E-06;
    333. 

  333.       Rbc(5,10)   = -3.690408440000000E-06; Zbs(5,10)   = -4.748117770000000E-06;
    334. 

  334.       Rbc(6,10)   =  4.171438220000000E-06; Zbs(6,10)   =  7.838792440000000E-07;
    335. 

  335.       Rbc(7,10)   =  1.512182800000000E-06; Zbs(7,10)   =  1.488982220000000E-06;
    336. 

  336.       Rbc(8,10)   = -9.929365099999999E-07; Zbs(8,10)   = -3.261759890000000E-06;
    337. 

  337.       Rbc(9,10)   = -1.627494720000000E-06; Zbs(9,10)   = -4.826840650000000E-06;
    338. 

  338.       Rbc(10,10)  =  2.288831480000000E-06; Zbs(10,10)  = -2.453483540000000E-06;
    339. 

  339.     }
    340. 

  340.     X0_.ReInit(COORD_DIM * Nt0 * Np0);
    341. 

  341.     for (long t = 0; t < Nt0; t++) {
    342. 

  342.       for (long p = 0; p < Np0; p++) {
    343. 

  343.         Real theta = t * 2 * sctl::const_pi<Real>() / Nt0;
    344. 

  344.         Real phi   = p * 2 * sctl::const_pi<Real>() / Np0;
    345. 

  345.         Real R = 0, Z = 0;
    346. 

  346.         for (long i = -10; i <= 10; i++) {
    347. 

  347.           for (long j = -10; j <= 10; j++) {
    348. 

  348.             R += Rbc(i,j) * sctl::cos(j*phi - 5*i*theta);
    349. 

  349.             Z += Zbs(i,j) * sctl::sin(j*phi - 5*i*theta);
    350. 

  350.           }
    351. 

  351.         }
    352. 

  352.         X0_[0*Nt0*Np0 + t*Np0 + p] = R * sctl::cos(theta);
    353. 

  353.         X0_[1*Nt0*Np0 + t*Np0 + p] = R * sctl::sin(theta);
    354. 

  354.         X0_[2*Nt0*Np0 + t*Np0 + p] = Z;
    355. 

  355.       }
    356. 

  356.     }
    357. 

  357.   } else {
    358. 

  358.     SCTL_ASSERT_MSG(false, "Unknown surface type");
    359. 

  359.   }
    360. 

  360. }
    361. 

  361. 

  362. }
    363.