source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2010/Yukawa_SQ_v40.ipf @ 770

Last change on this file since 770 was 770, checked in by srkline, 12 years ago

Documentation additions and help links to be rolled into the 7.04 release. There still is a little more documentaiton to finish in the ModelDocs?.

File size: 88.8 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4/////////////////////////////////////////////
5//
6// One-Yukawa and Two-Yukawa strucutre factors
7//      Yun Liu, Wei-Ren Chen, and Sow-Hsin Chen, J. Chem. Phys. 122 (2005) 044507.
8//
9//
10// Converted from Matlab to C by Marcus Hennig on 5/12/10
11//
12// Converted to Igor XOP - SRK July 2010
13// -- There are many external calls and allocation/deallocation of memory, so the XOP is NOT THREADED
14// -- The function calculation is inherently AAO, so this XOP definition is DIFFERENT than
15//              all of the standard fitting functions.
16// -- so be sure that the P*S implementations are not threaded - although P(q) can be threaded
17//
18// *** passing in Z values of zero can cause the XOP to crash. test for them here and send good values.
19// -- the XOP will be modified to handle this and noted here when it is done. 0.001 seems to be OK
20//    as a low value.
21// -- for OneYukawa, 0.1 seems to be a reasonable minimum
22//
23// - remember that the dimensionless Q variable is Q*diameter
24//
25//
26// conversion to Igor from the c-code was not terribly painful, and very useful for debugging.
27//
28//
29//
30// as of September 2010:
31//
32// the one-component has not been tested at all
33//
34// -- the two component result nearly matches the result that Yun gets. I do need to relax the criteria for
35// rejecting solutions, however. The XOP code rejects solutions that Yun considers "good". I guess I
36// need all of the intermediate values (polynomial coefficients, solution vectors, etc.). Other than some of the
37// numerical values not matching up - the output S(q) looks to be correct.
38//
39// -- also, for some cases, the results are VERY finicky - ususally there is a threshold value say, in Z, where
40// going beyond that value is unstable. Here, in can be a bit random as to which values works and which do not.
41// It must be hitting some strange zeros in the functions.
42//
43//
44//              TO ADD:
45//
46// - a mechanism for plotting the potential, so that users have a good handle on what the parameters actually mean.
47// - as with plotting the potential, transform the K and Z into some energy scale if possible for interpretation.
48//
49//
50//
51//
52//
53/////////////////////////////////////////////
54
55
56
57
58Proc PlotOneYukawa(num,qmin,qmax)
59        Variable num=200,qmin=0.001,qmax=0.5
60        Prompt num "Enter number of data points for model: "
61        Prompt qmin "Enter minimum q-value (A^-1) for model: "
62        Prompt qmax "Enter maximum q-value (A^-1) for model: "
63       
64        Make/O/D/n=(num) xwave_1yuk,ywave_1yuk
65        xwave_1yuk = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))   
66        Make/O/D coef_1yuk = {0.1,50,-1,10}
67        make/o/t parameters_1yuk = {"volume fraction","Radius (A)","scale, K","charge, Z"}
68        Edit parameters_1yuk,coef_1yuk
69        Variable/G root:g_1yuk
70        g_1yuk := OneYukawa(coef_1yuk,ywave_1yuk,xwave_1yuk)
71//      g_1yuk := OneYukawaX(coef_1yuk,xwave_1yuk,ywave_1yuk)           //be sure to have x and y in the correct order
72        Display ywave_1yuk vs xwave_1yuk
73        ModifyGraph marker=29,msize=2,mode=4
74        Label bottom "q (A\\S-1\\M)"
75        Label left "Structure Factor"
76        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
77       
78        AddModelToStrings("OneYukawa","coef_1yuk","parameters_1yuk","1yuk")
79End
80
81//AAO version
82Function OneYukawa(cw,yw,xw) : FitFunc
83        Wave cw,yw,xw
84
85        if(abs(cw[3]) < 0.1)
86                cw[3] = 0.1
87        endif   
88       
89#if exists("OneYukawaX")               
90        OneYukawaX(cw,xw,yw)
91#else
92        yw = 0
93#endif
94        return(0)
95End
96
97
98// no igor code, return 0
99//
100Function fOneYukawa(w,x) : FitFunc
101        Wave w
102        Variable x
103                       
104        return (0)
105End
106
107//////////////////////////////////////////////////////////////
108Proc PlotTwoYukawa(num,qmin,qmax)
109        Variable num=200,qmin=0.001,qmax=0.5
110        Prompt num "Enter number of data points for model: "
111        Prompt qmin "Enter minimum q-value (A^-1) for model: "
112        Prompt qmax "Enter maximum q-value (A^-1) for model: "
113       
114        declare2YGlobals()              //only necessary if Igor code is used. Not needed if XOP code is used.
115       
116        Make/O/D/n=(num) xwave_2yuk,ywave_2yuk
117        xwave_2yuk = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))   
118        Make/O/D coef_2yuk = {0.2,50,6,10,-1,2}
119        make/o/t parameters_2yuk = {"volume fraction","Radius (A)","scale, K1","charge, Z1","scale, K2","charge, Z2"}
120        Edit parameters_2yuk,coef_2yuk
121        Variable/G root:g_2yuk
122        g_2yuk := TwoYukawa(coef_2yuk,ywave_2yuk,xwave_2yuk)
123        Display ywave_2yuk vs xwave_2yuk
124        ModifyGraph marker=29,msize=2,mode=4
125        Label bottom "q (A\\S-1\\M)"
126        Label left "Structure Factor"
127        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
128       
129        AddModelToStrings("TwoYukawa","coef_2yuk","parameters_2yuk","2yuk")
130       
131End
132
133
134//AAO version
135Function TwoYukawa(cw,yw,xw) : FitFunc
136        Wave cw,yw,xw
137
138        if(abs(cw[2]) < 0.001)
139                cw[2] = 0.001
140        endif
141        if(abs(cw[3]) < 0.001)
142                cw[3] = 0.001
143        endif
144        if(abs(cw[4]) < 0.001)
145                cw[4] = 0.001
146        endif
147        if(abs(cw[5]) < 0.001)
148                cw[5] = 0.001
149        endif   
150       
151       
152#if exists("TwoYukawaX")
153        TwoYukawaX(cw,xw,yw)
154#else
155        fTwoYukawa(cw,xw,yw)
156#endif
157        return(0)
158End
159
160Proc TestTheIgor2YUK()
161        //if the regular 2-yukawa procedure is already plotted
162        // -- then append it to thte graph yourself
163        Duplicate/O ywave_2yuk ywave_2yuk_Igor
164        Variable/G root:g_2yuk_Igor=0
165        g_2yuk_Igor := fTwoYukawa(coef_2yuk,xwave_2yuk,ywave_2yuk_Igor)
166End
167
168//
169Function fTwoYukawa(cw,xw,yw) : FitFunc
170        Wave cw,xw,yw
171
172        Variable Z1, Z2, K1, K2, phi,radius
173        phi = cw[0]
174        radius = cw[1]
175        K1 = cw[2]
176        Z1 = cw[3]
177        K2 = cw[4]
178        Z2 = cw[5]
179       
180        Variable a,b,c1,c2,d1,d2
181       
182        Variable ok,check,prnt
183        prnt = 0                //print out intermediates
184       
185        ok = TY_SolveEquations( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2, prnt )               // a,b,c1,c2,d1,d2 are returned
186        if(ok)
187                check = TY_CheckSolution( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
188                if(prnt)
189                        printf "solution = (%g, %g, %g, %g, %g, %g) check = %d\r", a, b, c1, c2, d1, d2, check
190                endif
191
192//              if(check)
193                if(ok)                          //if(ok) simply takes the best solution, not necessarily one that passes TY_CheckSolution
194                        yw = SqTwoYukawa(xw*radius*2, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2)
195//                      printf("%g      %g\n",q,sq)
196                endif
197        endif
198     
199        return (0)
200End
201
202
203///////////////////// converted procedures from c-code ////////////////////////////
204
205/// there were two functions defined as TY_q: one as TY_Q and one as TY_q. I renamed the TY_Q function as TY_capQ, and left TY_q unchanged
206
207// function TY_W change to TY_capW, since there is a wave named TY_w
208
209
210
211
212
213Static Function chop(x)
214        Variable x
215
216        if ( abs(x) < 1e-6 )
217                return 0
218        else
219                return x
220        endif
221       
222end
223
224Static Function pow(a,b)
225        Variable a,b
226       
227        return (a^b)
228end
229
230///*
231// ==================================================================================================
232//
233// The two-yukawa structure factor is uniquley determined by 6 parameters a, b, c1, c2, d1, d2,
234// which are the solution of a system of 6 equations ( 4 linear, 2 nonlinear ). The solution can
235// constructed by the roots of a polynomial of 22nd degree. For more details see attached
236// Mathematica notebook, where a derivation is given
237//
238// ==================================================================================================
239// */
240
241// these all may need to be declared as global variables !!
242//
243// - they are defined in a global scope in the c-code!
244//
245// - change the data folder
246Function declare2YGlobals()
247
248        NewDataFolder/O/S root:yuk
249       
250        Variable/G TY_q22
251        Variable/G TY_qa12, TY_qa21, TY_qa22, TY_qa23, TY_qa32
252        Variable/G TY_qb12, TY_qb21, TY_qb22, TY_qb23, TY_qb32
253        Variable/G TY_qc112, TY_qc121, TY_qc122, TY_qc123, TY_qc132
254        Variable/G TY_qc212, TY_qc221, TY_qc222, TY_qc223, TY_qc232
255        Variable/G TY_A12, TY_A21, TY_A22, TY_A23, TY_A32, TY_A41, TY_A42, TY_A43, TY_A52
256        Variable/G TY_B12, TY_B14, TY_B21, TY_B22, TY_B23, TY_B24, TY_B25, TY_B32, TY_B34
257        Variable/G TY_F14, TY_F16, TY_F18, TY_F23, TY_F24, TY_F25, TY_F26, TY_F27, TY_F28, TY_F29, TY_F32, TY_F33, TY_F34, TY_F35, TY_F36, TY_F37, TY_F38, TY_F39, TY_F310
258        Variable/G TY_G13, TY_G14, TY_G15, TY_G16, TY_G17, TY_G18, TY_G19, TY_G110, TY_G111, TY_G112, TY_G113, TY_G22, TY_G23, TY_G24, TY_G25, TY_G26, TY_G27, TY_G28, TY_G29, TY_G210, TY_G211, TY_G212, TY_G213, TY_G214
259       
260        SetDataFolder root:
261        //this is an array, already global TY_w[23];
262
263End
264
265
266Function TY_sigma( s, Z1,  Z2, a,  b,  c1,  c2,  d1,  d2 )
267        Variable  s, Z1,  Z2, a,  b,  c1,  c2,  d1,  d2
268
269        return -(a / 2. + b + c1 * exp( -Z1 ) + c2 * exp( -Z2 )) / s + a * pow( s, -3 ) + b * pow( s, -2 ) + ( c1 + d1 ) * pow( s + Z1, -1 ) + ( c2 + d2 ) * pow( s + Z2, -1 )
270end
271
272Function TY_tau(  s, Z1,  Z2, a,  b,  c1,  c2 )
273        Variable   s, Z1,  Z2, a,  b,  c1,  c2
274       
275        return b * pow( s, -2 ) + a * ( pow( s, -3 ) + pow( s, -2 ) ) - pow( s, -1 ) * ( c1 * Z1 * exp( -Z1 ) * pow( s + Z1, -1 ) + c2 * Z2 * exp( -Z2 ) * pow( s + Z2, -1 ) )
276end
277
278Function TY_q(  s, Z1,  Z2, a,  b,  c1,  c2,  d1,  d2 )
279        Variable  s, Z1,  Z2, a,  b,  c1,  c2,  d1,  d2
280        return TY_sigma(s, Z1, Z2, a, b, c1, c2, d1, d2) - exp( -s ) * TY_tau(s, Z1, Z2, a,b, c1, c2)
281end
282
283Function TY_g(  s, phi,  Z1,  Z2, a,  b,  c1,  c2,  d1,  d2 )
284        Variable   s, phi,  Z1,  Z2, a,  b,  c1,  c2,  d1,  d2
285        return s * TY_tau( s, Z1, Z2, a, b, c1, c2 ) * exp( -s ) / ( 1 - 12 * phi * TY_q( s, Z1, Z2, a, b, c1, c2, d1, d2 ) )
286end
287
288///*
289// ==================================================================================================
290//
291// Structure factor for the potential
292//
293// V(r) = -kB * T * ( K1 * exp[ -Z1 * (r - 1)] / r + K2 * exp[ -Z2 * (r - 1)] / r ) for r > 1
294// V(r) = inf for r <= 1
295//
296// The structure factor is parametrized by (a, b, c1, c2, d1, d2)
297// which depend on (K1, K2, Z1, Z2, phi). 
298//
299// ==================================================================================================
300// */
301
302Function TY_hq(  q,  Z,  K,  v )
303        Variable   q,  Z,  K,  v
304       
305        if ( q == 0)
306                return (exp(-2.*Z)*(v + (v*(-1. + Z) - 2.*K*Z)*exp(Z))*(-(v*(1. + Z)) + (v + 2.*K*Z*(1. + Z))*exp(Z))*pow(K,-1)*pow(Z,-4))/4.
307        else
308       
309                variable t1, t2, t3, t4
310               
311                t1 = ( 1. - v / ( 2. * K * Z * exp( Z ) ) ) * ( ( 1. - cos( q ) ) / ( q*q ) - 1. / ( Z*Z + q*q ) )
312                t2 = ( v*v * ( q * cos( q ) - Z * sin( q ) ) ) / ( 4. * K * Z*Z * q * ( Z*Z + q*q ) )
313                t3 = ( q * cos( q ) + Z * sin( q ) ) / ( q * ( Z*Z + q*q ) )
314                t4 = v / ( Z * exp( Z ) ) - v*v / ( 4. * K * Z*Z * exp( 2. * Z ) ) - K
315               
316                return v / Z * t1 - t2 + t3 * t4
317        endif
318end
319
320
321Function TY_pc(  q, Z1,  Z2, K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2 )
322        Variable   q, Z1,  Z2, K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2
323       
324        variable v1 = 24. * phi * K1 * exp( Z1 ) * TY_g( Z1, phi, Z1, Z2, a, b, c1, c2, d1, d2 )
325        variable v2 = 24. * phi * K2 * exp( Z2 ) * TY_g( Z2, phi, Z1, Z2, a, b, c1, c2, d1, d2 )
326       
327        variable a0 = a * a
328        variable b0 = -12. * phi *( pow( a + b,2 ) / 2. + a * ( c1 * exp( -Z1 ) + c2 * exp( -Z2 ) ) )
329       
330        variable t1, t2, t3
331       
332        if ( q == 0 )
333                t1 = a0 / 3.
334                t2 = b0 / 4.
335                t3 = a0 * phi / 12.
336        else
337                t1 = a0 * ( sin( q ) - q * cos( q ) ) / pow( q, 3 )
338                t2 = b0 * ( 2. * q * sin( q ) - ( q * q - 2. ) * cos( q ) - 2. ) / pow( q, 4 )
339                t3 = a0 * phi * ( ( q*q - 6. ) * 4. * q * sin( q ) - ( pow( q, 4 ) - 12. * q*q + 24.) * cos( q ) + 24. ) / ( 2. * pow( q, 6 ) )
340        endif
341       
342        variable t4 = TY_hq( q, Z1, K1, v1 ) + TY_hq( q, Z2, K2, v2 )
343       
344        return -24. * phi * ( t1 + t2 + t3 + t4 )
345end
346
347Function SqTwoYukawa(  q, Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2 )
348        variable   q, Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2
349       
350        if ( Z1 == Z2 )
351                // one-yukawa potential
352                return 0
353        else
354                // two-yukawa potential
355                return 1. / ( 1. - TY_pc( q, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
356        endif
357end
358
359///*
360//==================================================================================================
361//
362// Non-linear eqaution system that determines the parameter for structure factor
363// 
364//==================================================================================================
365//*/
366
367Function TY_LinearEquation_1(  Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2 )
368        Variable   Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2
369       
370        return b - 12. * phi * ( -a / 8. - b / 6. + d1 * pow( Z1, -2 ) + c1 * ( pow( Z1, -2 )  - exp( -Z1 ) * ( 0.5 + ( 1. + Z1 ) * pow( Z1, -2 ) ) ) + d2 * pow( Z2, -2 ) + c2 * ( pow( Z2, -2 ) - exp( -Z2 )* ( 0.5 + ( 1. + Z2 ) * pow( Z2, -2 ) ) ) )
371end
372
373Function TY_LinearEquation_2(  Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2 )
374        Variable  Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2
375       
376        return 1. - a - 12. * phi * ( -a / 3. - b / 2. + d1 * pow( Z1, -1 ) + c1 * ( pow( Z1, -1 ) - ( 1. + Z1 ) * exp( -Z1 ) * pow( Z1, -1 ) ) + d2 * pow( Z2, -1 ) + c2 * ( pow( Z2, -1 ) - ( 1. + Z2 ) * exp( -Z2 ) * pow( Z2, -1 ) ) )
377end     
378
379Function TY_LinearEquation_3(  Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2 )
380        Variable  Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2
381                                                       
382        return K1 * exp( Z1 ) - d1 * Z1 * ( 1. - 12. * phi * TY_q( Z1, Z1, Z2, a, b, c1, c2, d1, d2 ) )
383end
384
385Function TY_LinearEquation_4(  Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2 )
386        Variable  Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2
387       
388        return K2 * exp( Z2 ) - d2 * Z2 * ( 1. - 12. * phi * TY_q( Z2, Z1, Z2, a, b, c1, c2, d1, d2 ) )
389end
390
391Function TY_NonlinearEquation_1(  Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2 )
392        Variable  Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2
393       
394        return c1 + d1 - 12. * phi * ( ( c1 + d1 ) * TY_sigma( Z1, Z1, Z2, a, b, c1, c2, d1, d2 ) - c1 * TY_tau( Z1, Z1, Z2, a, b, c1, c2 ) * exp( -Z1 ) )
395end
396
397Function TY_NonlinearEquation_2(  Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2 )
398        Variable  Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2
399       
400        return c2 + d2 - 12. * phi * ( ( c2 + d2 ) * TY_sigma( Z2, Z1, Z2, a, b, c1, c2, d1, d2 ) - c2 * TY_tau( Z2, Z1, Z2, a, b, c1, c2 ) * exp( -Z2 ) )
401end
402
403// Check the computed solutions satisfy the system of equations
404Function TY_CheckSolution(  Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2 )
405        variable   Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2
406       
407        variable eq_1 = chop( TY_LinearEquation_1( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
408        variable eq_2 = chop( TY_LinearEquation_2( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
409        variable eq_3 = chop( TY_LinearEquation_3( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
410        variable eq_4 = chop( TY_LinearEquation_4( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
411        variable eq_5 = chop( TY_NonlinearEquation_1( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
412        variable eq_6 = chop( TY_NonlinearEquation_2( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
413       
414//      printf("Check of solution = %g %g %g %g %g %g\r",eq_1,eq_2,eq_3,eq_4,eq_5,eq_6);
415        // check if all equation are zero
416        return ( eq_1 == 0 && eq_2 == 0 && eq_3 == 0 && eq_4 == 0 && eq_5 == 0 && eq_6 == 0 )
417end
418
419Function TY_ReduceNonlinearSystem( Z1, Z2,  K1,  K2,  phi,  prnt )
420        Variable  Z1, Z2,  K1,  K2,  phi,  prnt
421       
422       
423//      /* solution of the 4 linear equations depending on d1 and d2, the solution is polynomial
424//       in d1, d2. We represend the solution as determiants obtained by Cramer's rule
425//       which can be expressed by their coefficient matrices
426//       */
427       
428        Variable m11 = (3.*phi)/2.
429        Variable m13 = 6.*phi*exp(-Z1)*(2. + Z1*(2. + Z1) - 2.*exp(Z1))*pow(Z1,-2)
430        Variable m14 = 6.*phi*exp(-Z2)*(2. + Z2*(2. + Z2) - 2.*exp(Z2))*pow(Z2,-2)
431        Variable m23 = -12.*phi*exp(-Z1)*(-1. - Z1 + exp(Z1))*pow(Z1,-1)
432        Variable m24 = -12.*phi*exp(-Z2)*(-1. - Z2 + exp(Z2))*pow(Z2,-1)
433        Variable m31 = -6.*phi*exp(-Z1)*pow(Z1,-2)*(2.*(1 + Z1) + exp(Z1)*(-2. + pow(Z1,2)))
434        Variable m32 = -12.*phi*(-1. + Z1 + exp(-Z1))*pow(Z1,-1)
435        Variable m33 = 6.*phi*exp(-2.*Z1)*pow(-1. + exp(Z1),2)
436        Variable m34 = 12.*phi*exp(-Z1 - Z2)*(Z2 - (Z1 + Z2)*exp(Z1) + Z1*exp(Z1 + Z2))*pow(Z1 + Z2,-1)
437        Variable m41 = -6.*phi*exp(-Z2)*pow(Z2,-2)*(2.*(1. + Z2) + exp(Z2)*(-2. + pow(Z2,2)))
438        Variable m42 = -12.*phi*(-1. + Z2 + exp(-Z2))*pow(Z2,-1)
439        Variable m43 = 12.*phi*exp(-Z1 - Z2)*(Z1 - (Z1 + Z2 - Z2*exp(Z1))*exp(Z2))*pow(Z1 + Z2,-1)
440        Variable m44 = 6.*phi*exp(-2*Z2)*pow(-1. + exp(Z2),2)
441       
442//      /* determinant of the linear system expressed as coefficient matrix in d1, d2 */
443       
444        NVAR TY_q22 = root:yuk:TY_q22
445       
446        TY_q22 = m14*(-(m33*m42) + m23*(m32*m41 - m31*m42) + m32*m43 + (4.*m11*(-3.*m33*m41 + 2.*m33*m42 + 3.*m31*m43 - 2.*m32*m43))/3.)
447        TY_q22 +=  m13*(m34*m42 + m24*(-(m32*m41) + m31*m42) - m32*m44 + (4.*m11*(3.*m34*m41 - 2.*m34*m42 - 3.*m31*m44 + 2.*m32*m44))/3.)
448        TY_q22 += (3.*m24*(m33*(3.*m41 + 4.*m11*m41 - 3.*m11*m42) + (-3.*m31 - 4.*m11*m31 + 3.*m11*m32)*m43) + 3.*m23*(-3.*m34*m41 - 4.*m11*m34*m41 + 3.*m11*m34*m42 + 3.*m31*m44 + 4.*m11*m31*m44 - 3.*m11*m32*m44) - (m34*m43 - m33*m44)*pow(3. - 2.*m11,2))/9.
449       
450        if( prnt )
451                printf "\rDet = \r"
452//              printf "%f\t%f\r%f\t%f\r", 0., 0., 0., TY_q22
453                printf "TY_q22 = %15.12g\r",TY_q22
454        endif
455       
456//      /* Matrix representation of the determinant of the of the system where row refering to
457//       the variable a is replaced by solution vector */
458       
459        NVAR TY_qa12 = root:yuk:TY_qa12
460        NVAR TY_qa21 = root:yuk:TY_qa21
461        NVAR TY_qa22 = root:yuk:TY_qa22
462        NVAR TY_qa23 = root:yuk:TY_qa23
463        NVAR TY_qa32 = root:yuk:TY_qa32
464
465        Variable t1,t2,t3,t4,t5,t6,t7,t8,t9,t10
466        Variable t11,t12,t13,t14,t15,t16,t17,t18,t19,t20                //simply to keep the line length small enough
467       
468        TY_qa12 = (K1*(3.*m14*(m23*m42 - 4.*m11*m43) - 3.*m13*(m24*m42 - 4.*m11*m44) + (3. + 4.*m11)*(m24*m43 - m23*m44))*exp(Z1))/3.
469       
470        TY_qa21 = -(K2*(3.*m14*(m23*m32 - 4.*m11*m33) - 3.*m13*(m24*m32 - 4.*m11*m34) + (3. + 4.*m11)*(m24*m33 - m23*m34))*exp(Z2))/3.
471       
472        TY_qa22 = m14*(-(m23*m42*Z1) + 4.*m11*m43*Z1 - m33*(m42 + 4.*m11*Z2) + m32*(m43 + m23*Z2)) + (3.*m13*(m24*m42*Z1 - 4.*m11*m44*Z1 + m34*(m42 + 4.*m11*Z2) - m32*(m44 + m24*Z2)) + (3. + 4.*m11)*(-(m24*m43*Z1) + m23*m44*Z1 - m34*(m43 + m23*Z2) + m33*(m44 + m24*Z2)))/3.
473       
474
475        t1 = (2.*(-3.*m13*m42 + 3.*m43 + 4.*m11*m43)*Z1*pow(Z2,2) - m33*(Z1 + Z2)*(6.*m42 + (3. + 4.*m11)*pow(Z2,2)) +  3.*m32*(Z1 + Z2)*(2.*m43 + m13*pow(Z2,2)))
476        t2 = (2.*(3.*m14*m42 - 3.*m44 - 4.*m11*m44)*Z1*pow(Z2,2) + m34*(Z1 + Z2)*(6.*m42 + (3. + 4.*m11)*pow(Z2,2)) - 3.*m32*(Z1 + Z2)*(2.*m44 + m14*pow(Z2,2)))
477        t3 = (3.*(m14*m33*m42 - m13*m34*m42 - m14*m32*m43 + m34*m43 + m13*m32*m44 - m33*m44)*Z2*(Z1 + Z2) +  2.*m11*(6.*(-(m14*m43) + m13*m44)*Z1*pow(Z2,2) + m34*(Z1 + Z2)*(2.*m43*(-3. + Z2) - 3.*m13*pow(Z2,2)) +  m33*(Z1 + Z2)*(6.*m44 - 2.*m44*Z2 + 3.*m14*pow(Z2,2))))
478                 
479        TY_qa23 = 2.*phi*pow(Z2,-2)*(m24*t1 + m23*t2 + 2.*t3)*pow(Z1 + Z2,-1)   
480       
481       
482       
483        t1 = ((-3.*m13*m42 + (3. + 4.*m11)*m43)*(Z1 + Z2)*pow(Z1,2) - 2.*m33*(3.*m42*(Z1 + Z2) + (3. + 4.*m11)*Z2*pow(Z1,2)) + 6.*m32*(m43*(Z1 + Z2) + m13*Z2*pow(Z1,2)))
484        t2 = ((3.*m14*m42 - (3. + 4.*m11)*m44)*(Z1 + Z2)*pow(Z1,2) + m34*(6.*m42*(Z1 + Z2) + 2.*(3. + 4.*m11)*Z2*pow(Z1,2)) - 6.*m32*(m44*(Z1 + Z2) + m14*Z2*pow(Z1,2)))
485        t3 = (3.*(m14*m33*m42 - m13*m34*m42 - m14*m32*m43 + m34*m43 + m13*m32*m44 - m33*m44)*Z1*(Z1 + Z2) + 2.*m11*(-3.*(m14*m43 - m13*m44)*(Z1 + Z2)*pow(Z1,2) + 2.*m34*(m43*(-3 + Z1)*(Z1 + Z2) - 3.*m13*Z2*pow(Z1,2)) + m33*(-2.*m44*(-3. + Z1)*(Z1 + Z2) + 6.*m14*Z2*pow(Z1,2))))
486       
487        TY_qa32 = 2.*phi*pow(Z1,-2)*(m24*t1 + m23*t2 + 2.*t3)*pow(Z1 + Z2,-1)
488               
489        if( prnt )
490                printf "\rDet_a = \r"
491//              printf  "%f\t%f\t%f\r%f\t%f\t%f\r%f\t%f\t%f\r", 0., TY_qa12, 0., TY_qa21, TY_qa22, TY_qa23,  0., TY_qa32, 0.
492                printf "TY_qa12 = %15.12g\r",TY_qa12
493                printf "TY_qa21 = %15.12g\r",TY_qa21
494                printf "TY_qa22 = %15.12g\r",TY_qa22
495                printf "TY_qa23 = %15.12g\r",TY_qa23
496                printf "TY_qa32 = %15.12g\r",TY_qa32
497        endif
498       
499//      /* Matrix representation of the determinant of the of the system where row refering to
500//       the variable b is replaced by solution vector */
501
502        NVAR TY_qb12 = root:yuk:TY_qb12
503        NVAR TY_qb21 = root:yuk:TY_qb21
504        NVAR TY_qb22 = root:yuk:TY_qb22
505        NVAR TY_qb23 = root:yuk:TY_qb23
506        NVAR TY_qb32 = root:yuk:TY_qb32
507
508        TY_qb12 = (K1*(-3.*m11*m24*m43 + m14*(-3.*m23*m41 + (-3. + 8.*m11)*m43) + 3.*m11*m23*m44 + m13*(3.*m24*m41 + 3.*m44 - 8.*m11*m44))*exp(Z1))/3.
509       
510        TY_qb21 = (K2*(-3.*m13*m24*m31 + 3.*m11*m24*m33 + m14*(3.*m23*m31 + (3. - 8.*m11)*m33) - 3.*m13*m34 + 8.*m11*m13*m34 - 3.*m11*m23*m34)*exp(Z2))/3.
511       
512        TY_qb22 = m13*(m31*m44 - m24*m41*Z1 - m44*Z1 + (8.*m11*m44*Z1)/3. + m24*m31*Z2 + m34*(-m41 + Z2 - (8.*m11*Z2)/3.)) + m14*(m23*m41*Z1 + m43*Z1 - (8.*m11*m43*Z1)/3. + m33*(m41 - Z2 + (8.*m11*Z2)/3.) - m31*(m43 + m23*Z2)) +  m11*(m24*m43*Z1 - m23*m44*Z1 + m34*(m43 + m23*Z2) - m33*(m44 + m24*Z2))   
513       
514        t1 = (-(m14*m33*m41) + m13*m34*m41 + m14*m31*m43 - m11*m34*m43 - m13*m31*m44 + m11*m33*m44)
515        t2 = (-3.*m11*m24*m43 + m14*(-3.*m23*m41 + (-3. + 8.*m11)*m43) + 3.*m11*m23*m44 + m13*(3.*m24*m41 + 3.*m44 - 8.*m11*m44))
516        t3 = (3.*m24*(m33*m41 - m31*m43) + m23*(-3.*m34*m41 + 3.*m31*m44) + (-3. + 8.*m11)*(m34*m43 - m33*m44))
517       
518        TY_qb23 = 2.*phi*(3.*m14*m23*m31 - 3.*m13*m24*m31 + 3.*m14*m33 - 8.*m11*m14*m33 + 3.*m11*m24*m33 - 3.*m13*m34 + 8.*m11*m13*m34 -  3.*m11*m23*m34 + 2.*t3*  pow(Z2,-2) + 6.*t1*pow(Z2,-1) +  2.*t2*Z1*pow(Z1 + Z2,-1))
519       
520       
521        t1 = (-(m34*(m23*m41 + m43)) + m24*(m33*m41 - m31*m43) + (m23*m31 + m33)*m44)
522        t2 = (-(m14*m33*m41) + m13*m34*m41 + m14*m31*m43 - m13*m31*m44)
523        t3 = (m14*(2.*m23*m31 + 2.*m33 - m23*m41 - m43) + m13*(-2.*m34 + m24*(-2.*m31 + m41) + m44))
524        t4 = (16.*m34*m43 - 16.*m33*m44 - 6.*m34*m43*Z1 + 6.*m33*m44*Z1 + (6.*m24*m33 - 3.*m24*m43 + 8.*m14*(-2.*m33 + m43) + (8.*m13 - 3.*m23)*(2.*m34 - m44))*pow(Z1,2))
525        t5 = (2.*m34*m43*(8. - 3.*Z1) + 2.*m33*m44*(-8. + 3.*Z1) + (8.*m14*m43 - 3.*m24*m43 - 8.*m13*m44 + 3.*m23*m44)*pow(Z1,2))
526       
527        TY_qb32 = 2.*phi*pow(Z1,-2)*(6.*t1 +  6.*t2*Z1 +  3.*t3*pow(Z1,2) + (m11*Z2*t4 + m11*Z1*t5)* pow(Z1 + Z2,-1) + 6.*(-(m14*(m23*m31 + m33)) + m13*(m24*m31 + m34))*pow(Z1,3)*pow(Z1 + Z2,-1))
528               
529               
530        if( prnt )
531                printf "\rDet_b = \r"
532//              printf "%f\t%f\t%f\r%f\t%f\t%f\r%f\t%f\t%f\r", 0., TY_qb12, 0., TY_qb21, TY_qb22, TY_qb23, 0., TY_qb32, 0.
533                printf "TY_qb12 = %15.12g\r",TY_qb12
534                printf "TY_qb21 = %15.12g\r",TY_qb21
535                printf "TY_qb22 = %15.12g\r",TY_qb22
536                printf "TY_qb23 = %15.12g\r",TY_qb23
537                printf "TY_qb32 = %15.12g\r",TY_qb32
538        endif
539       
540//      /* Matrix representation of the determinant of the of the system where row refering to
541//       the variable c1 is replaced by solution vector */
542        NVAR TY_qc112 = root:yuk:TY_qc112
543        NVAR TY_qc121 = root:yuk:TY_qc121
544        NVAR TY_qc122 = root:yuk:TY_qc122
545        NVAR TY_qc123 = root:yuk:TY_qc123
546        NVAR TY_qc132 = root:yuk:TY_qc132
547
548        TY_qc112 = -(K1*exp(Z1)*(9.*m24*m41 - 9.*m14*m42 + 3.*m11*(-12.*m14*m41 + 4.*m24*m41 + 8.*m14*m42 - 3.*m24*m42) + m44*pow(3. - 2.*m11,2)))/9.
549       
550        TY_qc121 = (K2*exp(Z2)*(9.*m24*m31 - 9.*m14*m32 + 3.*m11*(-12.*m14*m31 + 4.*m24*m31 + 8.*m14*m32 - 3.*m24*m32) + m34*pow(3. - 2.*m11,2)))/9.
551       
552        TY_qc122 = m14*(-4.*m11*m41*Z1 - m42*Z1 + (8.*m11*m42*Z1)/3. + m32*(-m41 + Z2 - (8.*m11*Z2)/3.) + m31*(m42 + 4.*m11*Z2)) + (3.*m34*((3. + 4.*m11)*m41 - 3.*m11*m42) + 9.*m11*m32*m44 + 9.*m24*m41*Z1 + 12.*m11*m24*m41*Z1 - 9.*m11*m24*m42*Z1 + 9.*m44*Z1 - 12.*m11*m44*Z1 + 9.*m11*m24*m32*Z2 - 3.*(3. + 4.*m11)*m31*(m44 + m24*Z2) - m34*Z2*pow(3. - 2.*m11,2) + 4.*m44*Z1*pow(m11,2))/9.
553       
554       
555        t1 = (m34*(Z1 + Z2)*(2.*m42 + Z2*(-2.*m41 + Z2)) - m32*(Z1 + Z2)*(2.*m44 + m14*Z2*(-2.*m41 + Z2)) - 2.*(m14*m42 - m44)*Z2*(-(Z1*Z2) + m31*(Z1 + Z2)))
556        t2 = (2.*(3.*m41 + 4.*m11*m41 - 3.*m11*m42)*Z1*pow(Z2,2) + 3.*m32*(Z1 + Z2)*(2.*m41 + m11*pow(Z2,2)) - m31*(Z1 + Z2)*(6.*m42 + (3. + 4.*m11)*pow(Z2,2)))
557        t3 = (8.*m42 + 4.*m41*(-3. + Z2) - 3.*m42*Z2 + 2.*pow(Z2,2))
558        t4 = (6.*m44 - 2.*m44*Z2 + 3.*m14*pow(Z2,2))
559        t5 = (-8.*m32*m44*Z1 + m32*m44*(-8. + 3.*Z1)*Z2 + (3.*m32*m44 - 4.*(m14*(m32 + 3.*m41 - 2.*m42) + m44)*Z1)*pow(Z2,2) +  m34*(Z1 + Z2)*t3 + 2.*m31*(Z1 + Z2)*t4 - 4.*m14*m32*pow(Z2,3))
560                       
561        TY_qc123 = (2.*phi*pow(Z2,-2)*(9.*t1 + 4.*(-2.*m44*Z1 + m34*(Z1 + Z2))*pow(m11,2)*pow(Z2,2) - 3.*m24*t2 - 6.*m11*t5)*pow(Z1 + Z2,-1))/3.;
562       
563       
564        t1 = ((m14*m42 - m44)*(2.*m31 - Z1)*Z1*(Z1 + Z2) - 2.*m34*(m42*(Z1 + Z2) - Z1*(-(Z1*Z2) + m41*(Z1 + Z2))) + 2.*m32*(m44*(Z1 + Z2) - m14*Z1*(-(Z1*Z2) + m41*(Z1 + Z2))))
565        t2 = (((3. + 4.*m11)*m41 - 3.*m11*m42)*(Z1 + Z2)*pow(Z1,2) + 6.*m32*(m41*(Z1 + Z2) + m11*Z2*pow(Z1,2)) - 2.*m31*(3.*m42*(Z1 + Z2) + (3. + 4.*m11)*Z2*pow(Z1,2)))
566        t3 = (-8.*m32*m44 + m34*(m42*(8. - 3.*Z1) + 4.*m41*(-3. + Z1)) - 4.*m31*m44*(-3. + Z1) + 3.*m32*m44*Z1 - 2.*(3.*m14*m41 - 2.*m14*m42 + m44)*pow(Z1,2))
567        t4 = (4.*(3.*m31 - 2.*m32)*m44 + Z1*(-4.*m31*m44 + 3.*m32*m44 - 2.*(m14*(-6.*m31 + 4.*m32 + 3.*m41 - 2.*m42) + m44)*Z1) + m34*(m42*(8. - 3.*Z1) + 4.*m41*(-3. + Z1) + 4.*pow(Z1,2)))
568       
569        TY_qc132 = (-2.*phi*pow(Z1,-2)*(9.*t1 + 4.*(-2.*m34*Z2 + m44*(Z1 + Z2))*pow(m11,2)*pow(Z1,2) +  3.*m24*t2 + 6.*m11*(Z1*t3 + Z2*t4))*pow(Z1 + Z2,-1))/3.;
570               
571               
572        if( prnt )
573                printf "\rDet_c1 = \r"
574//              printf "%f\t%f\t%f\r%f\t%f\t%f\r%f\t%f\t%f\r", 0., TY_qc112, 0., TY_qc121, TY_qc122, TY_qc123, 0., TY_qc132, 0.
575                printf "TY_qc112 = %15.12g\r",TY_qc112
576                printf "TY_qc121 = %15.12g\r",TY_qc121
577                printf "TY_qc122 = %15.12g\r",TY_qc122
578                printf "TY_qc123 = %15.12g\r",TY_qc123
579                printf "TY_qc132 = %15.12g\r",TY_qc132
580        endif
581       
582//      /* Matrix representation of the determinant of the of the system where row refering to
583//       the variable c1 is replaced by solution vector */
584        NVAR TY_qc212 = root:yuk:TY_qc212
585        NVAR TY_qc221 = root:yuk:TY_qc221
586        NVAR TY_qc222 = root:yuk:TY_qc222
587        NVAR TY_qc223 = root:yuk:TY_qc223
588        NVAR TY_qc232 = root:yuk:TY_qc232
589
590        TY_qc212 = (K1*exp(Z1)*(9*m23*m41 - 9*m13*m42 + 3*m11*(-12*m13*m41 + 4*m23*m41 + 8*m13*m42 - 3*m23*m42) + m43*pow(3 - 2*m11,2)))/9.
591       
592        TY_qc221 = -(K2*exp(Z2)*(9*m23*m31 - 9*m13*m32 + 3*m11*(-12*m13*m31 + 4*m23*m31 + 8*m13*m32 - 3*m23*m32) + m33*pow(3 - 2*m11,2)))/9.
593       
594        TY_qc222 = m13*(4*m11*m41*Z1 + m42*Z1 - (8*m11*m42*Z1)/3. + m32*(m41 - Z2 + (8*m11*Z2)/3.) - m31*(m42 + 4*m11*Z2)) + (9*m31*m43 - 9*(m23*m41 + m43)*Z1 + 9*m23*m31*Z2 + 3*m11*((-4*m23*m41 + 3*m23*m42 + 4*m43)*Z1 + 4*m31*(m43 + m23*Z2) - 3*m32*(m43 + m23*Z2)) + m33*(-3*(3 + 4*m11)*m41 + 9*m11*m42 + Z2*pow(3 - 2*m11,2)) - 4*m43*Z1*pow(m11,2))/9.
595       
596                       
597        t1 = (-(m33*(Z1 + Z2)*(2*m42 + Z2*(-2*m41 + Z2))) + m32*(Z1 + Z2)*(2*m43 + m13*Z2*(-2*m41 + Z2)) + 2*(m13*m42 - m43)*Z2*(-(Z1*Z2) + m31*(Z1 + Z2)))
598        t2 = (2*(3*m41 + 4*m11*m41 - 3*m11*m42)*Z1*pow(Z2,2) + 3*m32*(Z1 + Z2)*(2*m41 + m11*pow(Z2,2)) - m31*(Z1 + Z2)*(6*m42 + (3 + 4*m11)*pow(Z2,2)))
599        t3 = (-8*m32*m43*Z1 + m32*m43*(-8 + 3*Z1)*Z2 + (3*m32*m43 - 4*(m13*(m32 + 3*m41 - 2*m42) + m43)*Z1)*pow(Z2,2) + m33*(Z1 + Z2)*(8*m42 + 4*m41*(-3 + Z2) - 3*m42*Z2 + 2*pow(Z2,2)) + 2*m31*(Z1 + Z2)*(6*m43 - 2*m43*Z2 + 3*m13*pow(Z2,2)) - 4*m13*m32*pow(Z2,3))
600       
601        TY_qc223 = (2*phi*pow(Z2,-2)*(9*t1 - 4*(-2*m43*Z1 + m33*(Z1 + Z2))*pow(m11,2)*pow(Z2,2) + 3*m23*t2 + 6*m11*t3)*pow(Z1 + Z2,-1))/3.
602       
603       
604        t1 = ((m13*m42 - m43)*(2*m31 - Z1)*Z1*(Z1 + Z2) - 2*m33*(m42*(Z1 + Z2) - Z1*(-(Z1*Z2) + m41*(Z1 + Z2))) +       2*m32*(m43*(Z1 + Z2) - m13*Z1*(-(Z1*Z2) + m41*(Z1 + Z2))))
605        t2 = (((3 + 4*m11)*m41 - 3*m11*m42)*(Z1 + Z2)*pow(Z1,2) + 6*m32*(m41*(Z1 + Z2) + m11*Z2*pow(Z1,2)) -    2*m31*(3*m42*(Z1 + Z2) + (3 + 4*m11)*Z2*pow(Z1,2)))
606        t3 = (-8*m32*m43 + m33*(m42*(8 - 3*Z1) + 4*m41*(-3 + Z1)) - 4*m31*m43*(-3 + Z1) + 3*m32*m43*Z1 - 2*(3*m13*m41 - 2*m13*m42 + m43)*pow(Z1,2))
607        t4 = (4*(3*m31 - 2*m32)*m43 + Z1*(-4*m31*m43 + 3*m32*m43 - 2*(m13*(-6*m31 + 4*m32 + 3*m41 - 2*m42) + m43)*Z1) + m33*(m42*(8 - 3*Z1) + 4*m41*(-3 + Z1) + 4*pow(Z1,2)))
608       
609        TY_qc232 = (2*phi*pow(Z1,-2)*(9*t1 + 4*(-2*m33*Z2 + m43*(Z1 + Z2))*pow(m11,2)*pow(Z1,2) + 3*m23*t2 + 6*m11*(Z1*t3 + Z2*t4))*pow(Z1 + Z2,-1))/3.
610               
611               
612        if( prnt )
613                printf "\rDet_c2 = \r"
614//              printf "%f\t%f\t%f\r%f\t%f\t%f\r%f\t%f\t%f\r",  0., TY_qc212, 0.,  TY_qc221, TY_qc222, TY_qc223,  0., TY_qc232, 0.
615                printf "TY_qc212 = %15.12g\r",TY_qc212
616                printf "TY_qc221 = %15.12g\r",TY_qc221
617                printf "TY_qc222 = %15.12g\r",TY_qc222
618                printf "TY_qc223 = %15.12g\r",TY_qc223
619                printf "TY_qc232 = %15.12g\r",TY_qc232
620        endif
621       
622//      /* coefficient matrices of nonlinear equation 1 */
623        NVAR TY_A12 = root:yuk:TY_A12
624        NVAR TY_A21 = root:yuk:TY_A21
625        NVAR TY_A22 = root:yuk:TY_A22
626        NVAR TY_A23 = root:yuk:TY_A23
627        NVAR TY_A32 = root:yuk:TY_A32
628        NVAR TY_A41 = root:yuk:TY_A41
629        NVAR TY_A42 = root:yuk:TY_A42
630        NVAR TY_A43 = root:yuk:TY_A43
631        NVAR TY_A52 = root:yuk:TY_A52
632       
633        t1 = (Z1*(2*TY_qb12*(-1 + Z1)*(Z1 + Z2) - Z1*(2*TY_qc212*Z1 + TY_qc112*(Z1 + Z2))) + TY_qa12*(Z1 + Z2)*(-2 + pow(Z1,2)))
634        t2 = (exp(2*Z1)*t1 - TY_qc112*(Z1 + Z2)*pow(Z1,2) + 2*(Z1 + Z2)*exp(Z1)*(TY_qa12 + (TY_qa12 + TY_qb12)*Z1 + TY_qc112*pow(Z1,2)))
635                 
636        TY_A12 = 6*phi*TY_qc112*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*TY_qc212*exp(Z1)*(-Z2 + (Z1 + Z2)*exp(Z1))*pow(Z1,2) + exp(Z2)*t2)*pow(Z1 + Z2,-1);
637       
638       
639        t1 = (2*Z1*(TY_qb21*TY_qc112*(-1 + Z1)*(Z1 + Z2) + TY_qb12*TY_qc121*(-1 + Z1)*(Z1 + Z2) -  Z1*(TY_qc121*TY_qc212*Z1 + TY_qc112*(TY_qc121 + TY_qc221)*Z1 + TY_qc112*TY_qc121*Z2)) + TY_qa21*TY_qc112*(Z1 + Z2)*(-2 + pow(Z1,2)) + TY_qa12*TY_qc121*(Z1 + Z2)*(-2 + pow(Z1,2)))
640        t2 = (TY_qb21*TY_qc112 + TY_qc121*(TY_qa12 + TY_qb12 + 2*TY_qc112*Z1))
641        t3 = (2*(TY_qa12*TY_qc121 + TY_qa21*TY_qc112*(1 + Z1) + Z1*t2)*(Z1 + Z2)*exp(Z1) + exp(2*Z1)*t1 - 2*TY_qc112*TY_qc121*(Z1 + Z2)*pow(Z1,2))
642                 
643        TY_A21 = 6*phi*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*(TY_qc121*TY_qc212 + TY_qc112*TY_qc221)*exp(Z1)*(-Z2 + (Z1 + Z2)*exp(Z1))*pow(Z1,2) +  exp(Z2)*t3)*pow(Z1 + Z2,-1);
644       
645       
646        t1 = (TY_qb22*TY_qc112 + TY_qc122*(TY_qa12 + TY_qb12 + 2*TY_qc112*Z1))
647        t2 = (2*Z1*(TY_qb22*TY_qc112*(-1 + Z1)*(Z1 + Z2) + TY_qb12*TY_qc122*(-1 + Z1)*(Z1 + Z2) - Z1*(TY_qc122*TY_qc212*Z1 + TY_qc112*(TY_qc122 + TY_qc222)*Z1 + TY_qc112*TY_qc122*Z2)) + TY_qa22*TY_qc112*(Z1 + Z2)*(-2 + pow(Z1,2)) + TY_qa12*TY_qc122*(Z1 + Z2)*(-2 + pow(Z1,2)))
648        t3 = (12*phi*(TY_qa12*TY_qc122 + TY_qa22*TY_qc112*(1 + Z1) + Z1*t1)*(Z1 + Z2)*exp(Z1) - 2*phi*TY_qc112*TY_qc122*(Z1 + Z2)*pow(Z1,2) + exp(2*Z1)*(6*phi*t2 + TY_q22*TY_qc112*(Z1 + Z2)*pow(Z1,3)))
649                 
650        TY_A22 = exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(12*phi*(TY_qc122*TY_qc212 + TY_qc112*TY_qc222)*exp(Z1)*(-Z2 + (Z1 + Z2)*exp(Z1))*pow(Z1,2) +  exp(Z2)*t3)*pow(Z1 + Z2,-1);
651       
652       
653        t1 = ((TY_q22*TY_qc112 + TY_qc123*(TY_qc112 + TY_qc212) + TY_qc112*TY_qc223)*Z1 + TY_qc112*TY_qc123*Z2)
654        t2 = (TY_qa12*TY_qc123 + TY_qa23*TY_qc112*(1 + Z1) + Z1*(TY_qb23*TY_qc112 + TY_qc123*(TY_qa12 + TY_qb12 + 2*TY_qc112*Z1)))
655        t3 = (2*Z1*(TY_qb23*TY_qc112*(-1 + Z1)*(Z1 + Z2) + TY_qb12*TY_qc123*(-1 + Z1)*(Z1 + Z2) - Z1*t1) + TY_qa23*TY_qc112*(Z1 + Z2)*(-2 + pow(Z1,2)) +  TY_qa12*TY_qc123*(Z1 + Z2)*(-2 + pow(Z1,2)))
656       
657        TY_A23 = 6*phi*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*(TY_qc123*TY_qc212 + TY_qc112*TY_qc223)*exp(Z1)*(-Z2 + (Z1 + Z2)*exp(Z1))*pow(Z1,2) + exp(Z2)*(2*t2*(Z1 + Z2)*exp(Z1) +  exp(2*Z1)*t3 - 2*TY_qc112*TY_qc123*(Z1 + Z2)*pow(Z1,2)))*pow(Z1 + Z2,-1);
658       
659       
660        t1 = (TY_qb32*TY_qc112 + (TY_qa23 + TY_qb23)*TY_qc121 + (TY_qa21 + TY_qb21)*TY_qc123 + (TY_qa12 + TY_qb12)*TY_qc132 + TY_q22*TY_qc112*Z1 +  2*(TY_qc121*TY_qc123 + TY_qc112*TY_qc132)*Z1 + TY_qc122*(TY_qa22 + TY_qb22 + TY_qc122*Z1))
661        t2 = (TY_qc132*TY_qc212 + TY_qc123*TY_qc221 + TY_qc122*TY_qc222 + TY_qc121*TY_qc223 + TY_qc112*TY_qc232)
662        t3 = ((TY_q22 + TY_qc132)*TY_qc212 + TY_qc123*TY_qc221 + TY_qc122*TY_qc222 + TY_qc121*TY_qc223 + TY_qc112*TY_qc232)
663        t4 = (2*TY_qc121*TY_qc123 + 2*TY_qc112*TY_qc132 + pow(TY_qc122,2))
664        t5 = (6*phi*(2*Z1*(TY_qb12*(-1 + Z1)*(Z1 + Z2) - Z1*((TY_qc112 + TY_qc121 + TY_qc212)*Z1 + TY_qc112*Z2)) +  TY_qa12*(Z1 + Z2)*(-2 + pow(Z1,2))) + TY_qc122*(Z1 + Z2)*pow(Z1,3))
665        t6 = (-2*(TY_qa22*TY_qc122 + TY_qa21*TY_qc123 + TY_qa12*TY_qc132) - 2*(TY_qb32*TY_qc112 + TY_qb23*TY_qc121 + TY_qb22*TY_qc122 + TY_qb21*TY_qc123 + TY_qb12*TY_qc132)*Z1 +  (2*TY_qb32*TY_qc112 + 2*TY_qb23*TY_qc121 + (TY_qa22 + 2*TY_qb22 - TY_qc122)*TY_qc122 + (TY_qa21 + 2*TY_qb21 - 2*TY_qc121)*TY_qc123 + (TY_qa12 + 2*TY_qb12 - 2*TY_qc112)*TY_qc132)*pow(Z1,2))
666        t7 = -2*TY_qa22*TY_qc122*Z1 - 2*TY_qa21*TY_qc123*Z1 - 2*TY_qa12*TY_qc132*Z1 + TY_qa32*TY_qc112*(Z1 + Z2)*(-2 + pow(Z1,2)) + TY_qa23*TY_qc121*(Z1 + Z2)*(-2 + pow(Z1,2)) - 2*TY_qb32*TY_qc112*pow(Z1,2) - 2*TY_qb23*TY_qc121*pow(Z1,2) - 2*TY_qb22*TY_qc122*pow(Z1,2) - 2*TY_qb21*TY_qc123*pow(Z1,2) - 2*TY_qb12*TY_qc132*pow(Z1,2)
667        t8 = Z2*t6 + 2*TY_qb32*TY_qc112*pow(Z1,3) + 2*TY_qb23*TY_qc121*pow(Z1,3) + TY_qa22*TY_qc122*pow(Z1,3) + 2*TY_qb22*TY_qc122*pow(Z1,3) + TY_qa21*TY_qc123*pow(Z1,3) + 2*TY_qb21*TY_qc123*pow(Z1,3) - 2*TY_qc121*TY_qc123*pow(Z1,3) + TY_qa12*TY_qc132*pow(Z1,3)
668        t9 = (t7 + t8 + 2*TY_qb12*TY_qc132*pow(Z1,3) - 2*TY_qc112*TY_qc132*pow(Z1,3) - 2*TY_qc132*TY_qc212*pow(Z1,3) - 2*TY_qc123*TY_qc221*pow(Z1,3) - 2*TY_qc122*TY_qc222*pow(Z1,3) - 2*TY_qc121*TY_qc223*pow(Z1,3) - 2*TY_qc112*TY_qc232*pow(Z1,3) - pow(TY_qc122,2)*pow(Z1,3))
669        t10 = (12*phi*(TY_qa23*TY_qc121 + TY_qa22*TY_qc122 + TY_qa21*TY_qc123 + TY_qa12*TY_qc132 + TY_qa32*TY_qc112*(1 + Z1) + Z1*t1)*(Z1 + Z2)*exp(Z1 + Z2) - 12*phi*t2*Z2*exp(Z1)*pow(Z1,2) + 12*phi*t3*(Z1 + Z2)*exp(2*Z1)*pow(Z1,2) - 6*phi*(Z1 + Z2)*exp(Z2)*t4*pow(Z1,2) + exp(2*Z1 + Z2)*(TY_q22*t5 + 6*phi*t9)) 
670                 
671        TY_A32 = exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*t10*pow(Z1 + Z2,-1);
672       
673       
674        t1 = ((-(TY_qc132*TY_qc221) - TY_qc121*TY_qc232)*Z2 + ((TY_q22 + TY_qc132)*TY_qc221 + TY_qc121*TY_qc232)*(Z1 + Z2)*exp(Z1))
675        t2 = (TY_qa21*TY_qc132 + TY_qa32*TY_qc121*(1 + Z1) + Z1*(TY_qb32*TY_qc121 + (TY_qa21 + TY_qb21)*TY_qc132 + TY_qc121*(TY_q22 + 2*TY_qc132)*Z1))
676        t3 = (-2*(TY_qa32*TY_qc121 + TY_qa21*(TY_q22 + TY_qc132)) - 2*(TY_qb32*TY_qc121 + TY_qb21*(TY_q22 + TY_qc132))*Z1 + (TY_q22*(TY_qa21 + 2*TY_qb21 - 2*TY_qc121) + TY_qc121*(TY_qa32 + 2*TY_qb32 - 2*TY_qc132) + (TY_qa21 + 2*TY_qb21)*TY_qc132)*pow(Z1,2))
677        t4 = (-2*(TY_qa32*TY_qc121 + TY_qa21*(TY_q22 + TY_qc132)) - 2*(TY_qb32*TY_qc121 + TY_qb21*(TY_q22 + TY_qc132))*Z1 + (TY_qa32*TY_qc121 + 2*TY_qb32*TY_qc121 + TY_qa21*TY_qc132 + 2*TY_qb21*TY_qc132 - 2*TY_qc121*TY_qc132 + TY_q22*(TY_qa21 + 2*TY_qb21 - 2*TY_qc121 - 2*TY_qc221) - 2*TY_qc132*TY_qc221 - 2*TY_qc121*TY_qc232)*pow(Z1,2))
678       
679        TY_A41 = 6*phi*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*exp(Z1)*t1*pow(Z1,2) + exp(Z2)*(2*t2*(Z1 + Z2)*exp(Z1) - 2*TY_qc121*TY_qc132*(Z1 + Z2)*pow(Z1,2) + exp(2*Z1)*(Z2*t3 + Z1*t4)))*pow(Z1 + Z2,-1);
680       
681       
682        t1 = (TY_qb32*TY_qc122 + (TY_qa22 + TY_qb22)*TY_qc132 + TY_qc122*(TY_q22 + 2*TY_qc132)*Z1)
683        t2 = (TY_qc132*TY_qc222 + TY_qc122*TY_qc232)
684        t3 = ((TY_q22 + TY_qc132)*TY_qc222 + TY_qc122*TY_qc232)
685        t4 = (2*Z1*(TY_qb32*TY_qc122*(-1 + Z1)*(Z1 + Z2) + TY_qb22*TY_qc132*(-1 + Z1)*(Z1 + Z2) - Z1*(TY_qc132*TY_qc222*Z1 + TY_qc122*(TY_qc132 + TY_qc232)*Z1 + TY_qc122*TY_qc132*Z2)) + TY_qa32*TY_qc122*(Z1 + Z2)*(-2 + pow(Z1,2)) + TY_qa22*TY_qc132*(Z1 + Z2)*(-2 + pow(Z1,2)))
686        t5 = (6*phi*t4 + (Z1 + Z2)*pow(TY_q22,2)*pow(Z1,3) + TY_q22*(6*phi*(2*Z1*(TY_qb22*(-1 + Z1)*(Z1 + Z2) - Z1*((TY_qc122 + TY_qc222)*Z1 + TY_qc122*Z2)) + TY_qa22*(Z1 + Z2)*(-2 + pow(Z1,2))) + TY_qc132*(Z1 + Z2)*pow(Z1,3)))
687        t6 = (12*phi*(TY_qa22*TY_qc132 + TY_qa32*TY_qc122*(1 + Z1) + Z1*t1)*(Z1 + Z2)*exp(Z1 + Z2) - 12*phi*t2*Z2*exp(Z1)*pow(Z1,2) + 12*phi*t3*(Z1 + Z2)*exp(2*Z1)*pow(Z1,2) - 12*phi*TY_qc122*TY_qc132*(Z1 + Z2)*exp(Z2)*pow(Z1,2) + exp(2*Z1 + Z2)*t5)
688               
689        TY_A42 = exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*t6*pow(Z1 + Z2,-1);
690       
691       
692        t1 = ((TY_qc132*TY_qc223 + TY_qc123*TY_qc232)*Z2 - ((TY_q22 + TY_qc132)*TY_qc223 + TY_qc123*TY_qc232)*(Z1 + Z2)*exp(Z1))
693        t2 = (TY_qa23*TY_qc132 + TY_qa32*TY_qc123*(1 + Z1) + Z1*(TY_qb32*TY_qc123 + (TY_qa23 + TY_qb23)*TY_qc132 + TY_qc123*(TY_q22 + 2*TY_qc132)*Z1))
694        t3 = (2*TY_qa32*TY_qc123 + 2*TY_qa23*(TY_q22 + TY_qc132) + 2*(TY_qb32*TY_qc123 + TY_qb23*(TY_q22 + TY_qc132))*Z1 - (TY_q22*(TY_qa23 + 2*TY_qb23 - 2*TY_qc123) + TY_qc123*(TY_qa32 + 2*TY_qb32 - 2*TY_qc132) + (TY_qa23 + 2*TY_qb23)*TY_qc132)*pow(Z1,2))
695        t4 = (2*TY_qa32*TY_qc123 + 2*TY_qa23*(TY_q22 + TY_qc132) + 2*(TY_qb32*TY_qc123 + TY_qb23*(TY_q22 + TY_qc132))*Z1 + (-(TY_qa32*TY_qc123) - (TY_qa23 + 2*TY_qb23)*TY_qc132 + TY_q22*(-TY_qa23 + 2*(-TY_qb23 + TY_qc123 + TY_qc132 + TY_qc223)) + 2*(-(TY_qb32*TY_qc123) + TY_qc132*(TY_qc123 + TY_qc223) + TY_qc123*TY_qc232) + 2*pow(TY_q22,2))*pow(Z1,2))
696       
697        TY_A43 = -6*phi*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*exp(Z1)*t1*pow(Z1,2) + exp(Z2)*(-2*t2*(Z1 + Z2)*exp(Z1) + 2*TY_qc123*TY_qc132*(Z1 + Z2)*pow(Z1,2) + exp(2*Z1)*(Z2*t3 + Z1*t4)))*pow(Z1 + Z2,-1);
698       
699       
700        t1 = (TY_qc132*Z2 - (TY_q22 + TY_qc132)*(Z1 + Z2)*exp(Z1))
701        t2 = (Z1*(-2*TY_qb32*(-1 + Z1)*(Z1 + Z2) + Z1*((TY_q22 + TY_qc132 + 2*TY_qc232)*Z1 + (TY_q22 + TY_qc132)*Z2)) -  TY_qa32*(Z1 + Z2)*(-2 + pow(Z1,2)))
702        t3 = ((TY_q22 + TY_qc132)*exp(2*Z1)*t2 + (Z1 + Z2)*pow(TY_qc132,2)*pow(Z1,2) - 2*TY_qc132*(Z1 + Z2)*exp(Z1)*(TY_qa32 + (TY_qa32 + TY_qb32)*Z1 + (TY_q22 + TY_qc132)*pow(Z1,2)))
703       
704        TY_A52 = -6*phi*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*TY_qc232*exp(Z1)*t1*pow(Z1,2) + exp(Z2)*t3)*pow(Z1 + Z2,-1);
705       
706       
707        // normalize A
708//      /*double norm_A = sqrt(pow(TY_A52,2)+pow(TY_A43,2)+pow(TY_A42, 2)+pow(TY_A41, 2)+pow(TY_A32, 2)+
709//                                               pow(TY_A23,2)+pow(TY_A22,2)+pow(TY_A21, 2)+pow(TY_A12, 2));
710//      TY_A12 /= norm_A;
711//      TY_A21 /= norm_A;
712//      TY_A22 /= norm_A;
713//      TY_A23 /= norm_A;
714//      TY_A32 /= norm_A;
715//      TY_A41 /= norm_A;
716//      TY_A42 /= norm_A;
717//      TY_A43 /= norm_A;
718//      TY_A52 /= norm_A;*/
719       
720        if( prnt )
721                printf "\rNonlinear equation 1 = \r"
722//              printf "%f\t\t%f\t\t%f\r", 0.,   TY_A12, 0.
723//              printf "%f\t\t%f\t\t%f\r", TY_A21, TY_A22, TY_A23
724//              printf "%f\t\t%f\t\t%f\r",  0.,  TY_A32, 0.
725//              printf "%f\t\t%f\t\t%f\r", TY_A41, TY_A42, TY_A43
726//              printf "%f\t\t%f\t\t%f\r", 0.,   TY_A52, 0.             
727                printf "TY_A12 = %15.12g\r",TY_A12
728                printf "TY_A21 = %15.12g\r",TY_A21
729                printf "TY_A22 = %15.12g\r",TY_A22
730                printf "TY_A23 = %15.12g\r",TY_A23
731                printf "TY_A32 = %15.12g\r",TY_A32
732                printf "TY_A41 = %15.12g\r",TY_A41
733                printf "TY_A42 = %15.12g\r",TY_A42
734                printf "TY_A43 = %15.12g\r",TY_A43
735                printf "TY_A52 = %15.12g\r",TY_A52
736        endif
737       
738//      /* coefficient matrices of nonlinear equation 2 */
739        NVAR TY_B12 = root:yuk:TY_B12
740        NVAR TY_B14 = root:yuk:TY_B14
741        NVAR TY_B21 = root:yuk:TY_B21
742        NVAR TY_B22 = root:yuk:TY_B22
743        NVAR TY_B23 = root:yuk:TY_B23
744        NVAR TY_B24 = root:yuk:TY_B24
745        NVAR TY_B25 = root:yuk:TY_B25
746        NVAR TY_B32 = root:yuk:TY_B32
747        NVAR TY_B34 = root:yuk:TY_B34
748       
749       
750       
751        t1 = (TY_qa12*TY_qc221 + TY_qa21*TY_qc212*(1 + Z2) + Z2*(TY_qb21*TY_qc212 + TY_qc221*(TY_qa12 + TY_qb12 + 2*TY_qc212*Z2)))
752        t2 = (-(TY_qc121*TY_qc212) - TY_qc112*TY_qc221)
753        t3 = (TY_qb21*TY_qc212*(-1 + Z2)*(Z1 + Z2) + TY_qb12*TY_qc221*(-1 + Z2)*(Z1 + Z2) - Z2*(TY_qc212*TY_qc221*Z1 + TY_qc112*TY_qc221*Z2 + TY_qc212*(TY_qc121 + TY_qc221)*Z2))
754        t4 = (exp(Z1)*(2*Z2*t3 + TY_qa21*TY_qc212*(Z1 + Z2)*(-2 + pow(Z2,2)) + TY_qa12*TY_qc221*(Z1 + Z2)*(-2 + pow(Z2,2))) + 2*(TY_qc121*TY_qc212 + TY_qc112*TY_qc221)*(Z1 + Z2)*pow(Z2,2))
755       
756        TY_B12 = 6*phi*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-2*TY_qc212*TY_qc221*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 2*exp(Z2)*((Z1 + Z2)*t1*exp(Z1) + t2*Z1*pow(Z2,2)) + exp(2*Z2)*t4)*pow(Z1 + Z2,-1);
757       
758       
759       
760       
761        t1 = ((Z1 + Z2)*(TY_qa12*TY_qc223 + TY_qa23*TY_qc212*(1 + Z2) + Z2*(TY_qb23*TY_qc212 + (TY_qa12 + TY_qb12)*TY_qc223 + TY_qc212*(TY_q22 + 2*TY_qc223)*Z2))*exp(Z1) + (-(TY_qc123*TY_qc212) - TY_qc112*TY_qc223)*Z1*pow(Z2,2))
762        t2 = (TY_qc123*TY_qc212 + TY_qc112*(TY_q22 + TY_qc223))
763        t3 = (TY_qa23*TY_qc212 + TY_qa12*(TY_q22 + TY_qc223))
764        t4 = (TY_qa23*TY_qc212 + TY_qa12*(TY_q22 + TY_qc223) + (TY_qb23*TY_qc212 + TY_qb12*(TY_q22 + TY_qc223))*Z1)
765        t5 = (-2*(TY_qb23*TY_qc212 + TY_qb12*(TY_q22 + TY_qc223)) + (TY_q22*(TY_qa12 + 2*TY_qb12 - 2*TY_qc212) + TY_qc212*(TY_qa23 + 2*TY_qb23 - 2*TY_qc223) +  (TY_qa12 + 2*TY_qb12)*TY_qc223)*Z1)
766        t6 = (TY_q22*(TY_qa12 + 2*TY_qb12 - 2*TY_qc112 - 2*TY_qc212) + TY_qc212*(TY_qa23 + 2*TY_qb23 - 2*TY_qc123 - 2*TY_qc223) + (TY_qa12 + 2*TY_qb12 - 2*TY_qc112)*TY_qc223)
767       
768        TY_B14 = 6*phi*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-2*TY_qc212*TY_qc223*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 2*exp(Z2)*t1 +  exp(2*Z2)*(2*t2*(Z1 + Z2)*pow(Z2,2) + exp(Z1)*(-2*t3*Z1 - 2*t4*Z2 + t5*pow(Z2,2) + t6*pow(Z2,3))))*pow(Z1 + Z2,-1);
769       
770       
771       
772        t1 = (TY_qc221*(Z1 + Z2)*exp(Z1)*pow(Z2,2))
773        t2 = (exp(Z1)*(Z2*(2*TY_qb21*(-1 + Z2)*(Z1 + Z2) - Z2*(2*TY_qc121*Z2 + TY_qc221*(Z1 + Z2))) + TY_qa21*(Z1 + Z2)*(-2 + pow(Z2,2))) +  2*TY_qc121*(Z1 + Z2)*pow(Z2,2))
774        t3 = (-(TY_qc121*Z1*pow(Z2,2)) + (Z1 + Z2)*exp(Z1)*(TY_qa21 + (TY_qa21 + TY_qb21)*Z2 + TY_qc221*pow(Z2,2)))
775       
776        TY_B21 = 6*phi*TY_qc221*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-t1 +  exp(2*Z2)*t2 + 2*exp(Z2)*t3)*pow(Z1 + Z2,-1);
777       
778       
779       
780        t1 = (TY_qb22*TY_qc221 + TY_qc222*(TY_qa21 + TY_qb21 + 2*TY_qc221*Z2))
781        t2 = ((Z1 + Z2)*(TY_qa21*TY_qc222 + TY_qa22*TY_qc221*(1 + Z2) + Z2*t1)*exp(Z1) + (-(TY_qc122*TY_qc221) - TY_qc121*TY_qc222)*Z1*pow(Z2,2))
782        t3 = (TY_qc122*TY_qc221 + TY_qc121*TY_qc222)
783        t4 = (TY_qb22*TY_qc221*(-1 + Z2)*(Z1 + Z2) + TY_qb21*TY_qc222*(-1 + Z2)*(Z1 + Z2) - Z2*(TY_qc221*TY_qc222*Z1 + TY_qc121*TY_qc222*Z2 + TY_qc221*(TY_qc122 + TY_qc222)*Z2))
784        t5 = (6*phi*(2*Z2*t4 + TY_qa22*TY_qc221*(Z1 + Z2)*(-2 + pow(Z2,2)) + TY_qa21*TY_qc222*(Z1 + Z2)*(-2 + pow(Z2,2))) + TY_q22*TY_qc221*(Z1 + Z2)*pow(Z2,3))
785       
786        TY_B22 = exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-12*phi*TY_qc221*TY_qc222*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 12*phi*exp(Z2)*t2 + exp(2*Z2)*(12*phi*t3*(Z1 + Z2)*pow(Z2,2) +  exp(Z1)*t5))*pow(Z1 + Z2,-1);
787       
788       
789       
790       
791        t1 = (2*TY_qc221*TY_qc223 + 2*TY_qc212*TY_qc232 + pow(TY_qc222,2))
792        t2 = (TY_qb32*TY_qc212 + (TY_qa23 + TY_qb23)*TY_qc221 + (TY_qa22 + TY_qb22)*TY_qc222 + (TY_qa21 + TY_qb21)*TY_qc223 + (TY_qa12 + TY_qb12)*TY_qc232 + Z2*(TY_q22*TY_qc221 + 2*TY_qc221*TY_qc223 + 2*TY_qc212*TY_qc232 + pow(TY_qc222,2)))
793        t3 = (-(TY_qc132*TY_qc212) - TY_qc123*TY_qc221 - TY_qc122*TY_qc222 - TY_qc121*TY_qc223 - TY_qc112*TY_qc232)
794        t4 = (TY_qc132*TY_qc212 + TY_qc123*TY_qc221 + TY_qc122*TY_qc222 + TY_qc121*(TY_q22 + TY_qc223) + TY_qc112*TY_qc232)
795        t5 = (TY_qa32*TY_qc212 + TY_qa23*TY_qc221 + TY_qa22*TY_qc222 + TY_qa21*(TY_q22 + TY_qc223) + TY_qa12*TY_qc232)
796        t6 = (TY_qa32*TY_qc212 + TY_qa23*TY_qc221 + TY_qa22*TY_qc222 + TY_qa21*(TY_q22 + TY_qc223) + TY_qa12*TY_qc232 + (TY_qb32*TY_qc212 + TY_qb23*TY_qc221 + TY_qb22*TY_qc222 + TY_qb21*(TY_q22 + TY_qc223) + TY_qb12*TY_qc232)*Z1)
797        t7 = (TY_qb32*TY_qc212 + TY_qb23*TY_qc221 + TY_qb22*TY_qc222 + TY_qb21*(TY_q22 + TY_qc223) + TY_qb12*TY_qc232)
798        t8 = (TY_q22*(TY_qa21 + 2*TY_qb21 - 2*TY_qc221) + (TY_qa22 + 2*TY_qb22 - TY_qc222)*TY_qc222 + TY_qc221*(TY_qa23 + 2*TY_qb23 - 2*TY_qc223) + TY_qa21*TY_qc223 + 2*TY_qb21*TY_qc223 + TY_qc212*(TY_qa32 + 2*TY_qb32 - 2*TY_qc232) + TY_qa12*TY_qc232 + 2*TY_qb12*TY_qc232)
799        t9 = (TY_qa21 + 2*TY_qb21 - 2*TY_qc121 - 2*TY_qc212 - 2*TY_qc221)
800        t10 = (TY_qa22 + 2*TY_qb22 - 2*TY_qc122 - TY_qc222)
801        t11 = (TY_qa23 + 2*TY_qb23 - 2*TY_qc123 - 2*TY_qc223)
802        t12 = (TY_qa32 + 2*TY_qb32 - 2*TY_qc132 - 2*TY_qc232)
803        t13 = ((Z1 + Z2)*exp(Z1)*(TY_qa23*TY_qc221 + TY_qa22*TY_qc222 + TY_qa21*TY_qc223 + TY_qa12*TY_qc232 + TY_qa32*TY_qc212*(1 + Z2) + Z2*t2) + t3*Z1*pow(Z2,2))
804        t14 = (TY_q22*t9 + t10*TY_qc222 + TY_qc221*t11 + (TY_qa21 + 2*TY_qb21 - 2*TY_qc121)*TY_qc223 + TY_qc212*t12 + (TY_qa12 + 2*TY_qb12 - 2*TY_qc112)*TY_qc232)
805        t15 = (-6*phi*(Z1 + Z2)*exp(Z1)*t1*pow(Z2,2) +  12*phi*exp(Z2)*t13 +  exp(2*Z2)*(12*phi*t4*(Z1 + Z2)*pow(Z2,2) +  exp(Z1)*(TY_q22*TY_qc222*(Z1 + Z2)*pow(Z2,3) - 6*phi*(2*t5*Z1 + 2*t6*Z2 - (-2*t7 +  t8*Z1)*pow(Z2,2) - t14*pow(Z2,3)))))
806       
807        TY_B23 = exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*t15*pow(Z1 + Z2,-1);
808       
809       
810        t1 = (TY_qa22*TY_qc223 + TY_qa23*TY_qc222*(1 + Z2) + Z2*(TY_qb23*TY_qc222 + (TY_qa22 + TY_qb22)*TY_qc223 + TY_qc222*(TY_q22 + 2*TY_qc223)*Z2))
811        t2 = (TY_qc123*TY_qc222 + TY_qc122*TY_qc223)
812        t3 = (TY_qc123*TY_qc222 + TY_qc122*(TY_q22 + TY_qc223))
813        t4 = (2*Z2*(TY_qb23*TY_qc222*(-1 + Z2)*(Z1 + Z2) + TY_qb22*TY_qc223*(-1 + Z2)*(Z1 + Z2) - Z2*(TY_qc222*TY_qc223*Z1 + TY_qc122*TY_qc223*Z2 + TY_qc222*(TY_qc123 + TY_qc223)*Z2)) + TY_qa23*TY_qc222*(Z1 + Z2)*(-2 + pow(Z2,2)) + TY_qa22*TY_qc223*(Z1 + Z2)*(-2 + pow(Z2,2)))
814        t5 = (6*phi*t4 + (Z1 + Z2)*pow(TY_q22,2)*pow(Z2,3) + TY_q22*(6*phi*(2*Z2*(TY_qb22*(-1 + Z2)*(Z1 + Z2) - Z2*(TY_qc222*Z1 + (TY_qc122 + TY_qc222)*Z2)) + TY_qa22*(Z1 + Z2)*(-2 + pow(Z2,2))) + TY_qc223*(Z1 + Z2)*pow(Z2,3)))
815        t6 = (12*phi*(Z1 + Z2)*t1*exp(Z1 + Z2) - 12*phi*TY_qc222*TY_qc223*(Z1 + Z2)*exp(Z1)*pow(Z2,2) - 12*phi*t2*Z1*exp(Z2)*pow(Z2,2) + 12*phi*t3*(Z1 + Z2)*exp(2*Z2)*pow(Z2,2) + exp(Z1 + 2*Z2)*t5)
816               
817        TY_B24 = exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*t6*pow(Z1 + Z2,-1);
818
819       
820        t1 = (exp(Z1)*(Z2*(-2*TY_qb23*(-1 + Z2)*(Z1 + Z2) + Z2*((TY_q22 + TY_qc223)*Z1 + (TY_q22 + 2*TY_qc123 + TY_qc223)*Z2)) -  TY_qa23*(Z1 + Z2)*(-2 + pow(Z2,2))) - 2*TY_qc123*(Z1 + Z2)*pow(Z2,2))
821        t2 = ((Z1 + Z2)*exp(Z1)*pow(TY_qc223,2)*pow(Z2,2) + (TY_q22 + TY_qc223)*exp(2*Z2)*t1 + 2*TY_qc223*exp(Z2)*(TY_qc123*Z1*pow(Z2,2) - (Z1 + Z2)*exp(Z1)*(TY_qa23 + (TY_qa23 + TY_qb23)*Z2 + (TY_q22 + TY_qc223)*pow(Z2,2))))
822       
823        TY_B25 = -6*phi*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*t2*pow(Z1 + Z2,-1);
824       
825       
826        t1 = (TY_qa21*TY_qc232 + TY_qa32*TY_qc221*(1 + Z2) + Z2*(TY_qb32*TY_qc221 + TY_qc232*(TY_qa21 + TY_qb21 + 2*TY_qc221*Z2)))
827        t2 = (-(TY_qc132*TY_qc221) - TY_qc121*TY_qc232)
828        t3 = (TY_qb32*TY_qc221*(-1 + Z2)*(Z1 + Z2) + TY_qb21*TY_qc232*(-1 + Z2)*(Z1 + Z2) - Z2*(TY_qc221*TY_qc232*Z1 + TY_qc121*TY_qc232*Z2 + TY_qc221*(TY_q22 + TY_qc132 + TY_qc232)*Z2))
829        t4 = (exp(Z1)*(2*Z2*t3 + TY_qa32*TY_qc221*(Z1 + Z2)*(-2 + pow(Z2,2)) + TY_qa21*TY_qc232*(Z1 + Z2)*(-2 + pow(Z2,2))) + 2*(TY_qc132*TY_qc221 + TY_qc121*TY_qc232)*(Z1 + Z2)*pow(Z2,2)) 
830       
831        TY_B32 = 6*phi*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-2*TY_qc221*TY_qc232*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 2*exp(Z2)*((Z1 + Z2)*t1*exp(Z1) + t2*Z1*pow(Z2,2)) + exp(2*Z2)*t4)*pow(Z1 + Z2,-1);
832       
833
834        t1 = (-((Z1 + Z2)*(TY_qa23*TY_qc232 + TY_qa32*TY_qc223*(1 + Z2) + Z2*(TY_qb32*TY_qc223 + TY_qc232*(TY_qa23 + TY_qb23 + TY_q22*Z2 + 2*TY_qc223*Z2)))*exp(Z1)) + (TY_qc132*TY_qc223 + TY_qc123*TY_qc232)*Z1*pow(Z2,2))
835        t2 = (TY_qc132*(TY_q22 + TY_qc223) + TY_qc123*TY_qc232)
836        t3 = (TY_qa32*(TY_q22 + TY_qc223) + TY_qa23*TY_qc232)
837        t4 = (TY_qa32*(TY_q22 + TY_qc223) + TY_qa23*TY_qc232 + (TY_qb32*(TY_q22 + TY_qc223) + TY_qb23*TY_qc232)*Z1)
838        t5 = (-2*(TY_qb32*(TY_q22 + TY_qc223) + TY_qb23*TY_qc232) + ((TY_qa32 + 2*TY_qb32)*(TY_q22 + TY_qc223) + (-2*TY_q22 + TY_qa23 + 2*TY_qb23 - 2*TY_qc223)*TY_qc232)*Z1)
839        t6 = (2*t3*Z1 + 2*t4*Z2 - t5*pow(Z2,2) + ((2*TY_q22 - TY_qa32 - 2*TY_qb32 + 2*TY_qc132)*(TY_q22 + TY_qc223) + (2*TY_q22 - TY_qa23 + 2*(-TY_qb23 + TY_qc123 + TY_qc223))*TY_qc232)*pow(Z2,3))
840                         
841        TY_B34 = -6*phi*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(2*TY_qc223*TY_qc232*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 2*exp(Z2)*t1 + exp(2*Z2)*(-2*t2*(Z1 + Z2)*pow(Z2,2) + exp(Z1)*t6))*pow(Z1 + Z2,-1);
842       
843
844//      /*double norm_B = sqrt(pow(TY_B12, 2)+pow(TY_B14, 2)+pow(TY_B21, 2)+pow(TY_B22, 2)+pow(TY_B23, 2)+pow(TY_B24, 2)+pow(TY_B25, 2)+pow(TY_B32, 2)+pow(TY_B34, 2));
845//     
846//      TY_B12 /= norm_B;
847//      TY_B14 /= norm_B;
848//      TY_B21 /= norm_B;
849//      TY_B22 /= norm_B;
850//      TY_B23 /= norm_B;
851//      TY_B24 /= norm_B;
852//      TY_B25 /= norm_B;
853//      TY_B32 /= norm_B;
854//      TY_B34 /= norm_B; */
855       
856        if( prnt )
857                printf "\rNonlinear equation 2 = \r"
858//              printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0.,  TY_B12, 0.,  TY_B14, 0.
859//              printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\r", TY_B21, TY_B22, TY_B23, TY_B24, TY_B25
860//              printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0.,  TY_B32, 0.,  TY_B34, 0.
861                printf "TY_B12 = %15.12g\r",TY_B12
862                printf "TY_B14 = %15.12g\r",TY_B14
863                printf "TY_B21 = %15.12g\r",TY_B21
864                printf "TY_B22 = %15.12g\r",TY_B22
865                printf "TY_B23 = %15.12g\r",TY_B23
866                printf "TY_B24 = %15.12g\r",TY_B24
867                printf "TY_B25 = %15.12g\r",TY_B25
868                printf "TY_B32 = %15.12g\r",TY_B32
869                printf "TY_B34 = %15.12g\r",TY_B34
870        endif
871       
872//      /* decrease order of nonlinear equation 1 by means of equation 2 */
873        NVAR TY_F14 = root:yuk:TY_F14
874        NVAR TY_F16 = root:yuk:TY_F16
875        NVAR TY_F18 = root:yuk:TY_F18
876        NVAR TY_F23 = root:yuk:TY_F23
877        NVAR TY_F24 = root:yuk:TY_F24
878        NVAR TY_F25 = root:yuk:TY_F25
879        NVAR TY_F26 = root:yuk:TY_F26
880        NVAR TY_F27 = root:yuk:TY_F27
881        NVAR TY_F28 = root:yuk:TY_F28
882        NVAR TY_F29 = root:yuk:TY_F29
883        NVAR TY_F32 = root:yuk:TY_F32
884        NVAR TY_F33 = root:yuk:TY_F33
885        NVAR TY_F34 = root:yuk:TY_F34
886        NVAR TY_F35 = root:yuk:TY_F35
887        NVAR TY_F36 = root:yuk:TY_F36
888        NVAR TY_F37 = root:yuk:TY_F37
889        NVAR TY_F38 = root:yuk:TY_F38
890        NVAR TY_F39 = root:yuk:TY_F39
891        NVAR TY_F310 = root:yuk:TY_F310
892       
893        TY_F14 = -(TY_A32*TY_B12*TY_B32) + TY_A52*pow(TY_B12,2) + TY_A12*pow(TY_B32,2)
894        TY_F16 = 2*TY_A52*TY_B12*TY_B14 - TY_A32*TY_B14*TY_B32 - TY_A32*TY_B12*TY_B34 + 2*TY_A12*TY_B32*TY_B34
895        TY_F18 = -(TY_A32*TY_B14*TY_B34) + TY_A52*pow(TY_B14,2) + TY_A12*pow(TY_B34,2)
896        TY_F23 = 2*TY_A52*TY_B12*TY_B21 - TY_A41*TY_B12*TY_B32 - TY_A32*TY_B21*TY_B32 + TY_A21*pow(TY_B32,2)
897        TY_F24 = 2*TY_A52*TY_B12*TY_B22 - TY_A42*TY_B12*TY_B32 - TY_A32*TY_B22*TY_B32 + TY_A22*pow(TY_B32,2)
898        TY_F25 = 2*TY_A52*TY_B14*TY_B21 + 2*TY_A52*TY_B12*TY_B23 - TY_A43*TY_B12*TY_B32 - TY_A41*TY_B14*TY_B32 - TY_A32*TY_B23*TY_B32 - TY_A41*TY_B12*TY_B34 - TY_A32*TY_B21*TY_B34 + 2*TY_A21*TY_B32*TY_B34 + TY_A23*pow(TY_B32,2)
899        TY_F26 = 2*TY_A52*TY_B14*TY_B22 + 2*TY_A52*TY_B12*TY_B24 - TY_A42*TY_B14*TY_B32 - TY_A32*TY_B24*TY_B32 - TY_A42*TY_B12*TY_B34 - TY_A32*TY_B22*TY_B34 + 2*TY_A22*TY_B32*TY_B34
900        TY_F27 = 2*TY_A52*TY_B14*TY_B23 + 2*TY_A52*TY_B12*TY_B25 - TY_A43*TY_B14*TY_B32 - TY_A32*TY_B25*TY_B32 - TY_A43*TY_B12*TY_B34 - TY_A41*TY_B14*TY_B34 - TY_A32*TY_B23*TY_B34 + 2*TY_A23*TY_B32*TY_B34 + TY_A21*pow(TY_B34,2)
901        TY_F28 = 2*TY_A52*TY_B14*TY_B24 - TY_A42*TY_B14*TY_B34 - TY_A32*TY_B24*TY_B34 + TY_A22*pow(TY_B34,2)
902        TY_F29 = 2*TY_A52*TY_B14*TY_B25 - TY_A43*TY_B14*TY_B34 - TY_A32*TY_B25*TY_B34 + TY_A23*pow(TY_B34,2)
903        TY_F32 = -(TY_A41*TY_B21*TY_B32) + TY_A52*pow(TY_B21,2)
904        TY_F33 = 2*TY_A52*TY_B21*TY_B22 - TY_A42*TY_B21*TY_B32 - TY_A41*TY_B22*TY_B32
905        TY_F34 = 2*TY_A52*TY_B21*TY_B23 - TY_A43*TY_B21*TY_B32 - TY_A42*TY_B22*TY_B32 - TY_A41*TY_B23*TY_B32 - TY_A41*TY_B21*TY_B34 + TY_A52*pow(TY_B22,2)
906        TY_F35 = 2*TY_A52*TY_B22*TY_B23 + 2*TY_A52*TY_B21*TY_B24 - TY_A43*TY_B22*TY_B32 - TY_A42*TY_B23*TY_B32 - TY_A41*TY_B24*TY_B32 - TY_A42*TY_B21*TY_B34 - TY_A41*TY_B22*TY_B34
907        TY_F36 = 2*TY_A52*TY_B22*TY_B24 + 2*TY_A52*TY_B21*TY_B25 - TY_A43*TY_B23*TY_B32 - TY_A42*TY_B24*TY_B32 - TY_A41*TY_B25*TY_B32 - TY_A43*TY_B21*TY_B34 - TY_A42*TY_B22*TY_B34 - TY_A41*TY_B23*TY_B34 + TY_A52*pow(TY_B23,2)
908        TY_F37 = 2*TY_A52*TY_B23*TY_B24 + 2*TY_A52*TY_B22*TY_B25 - TY_A43*TY_B24*TY_B32 - TY_A42*TY_B25*TY_B32 - TY_A43*TY_B22*TY_B34 - TY_A42*TY_B23*TY_B34 - TY_A41*TY_B24*TY_B34
909        TY_F38 = 2*TY_A52*TY_B23*TY_B25 - TY_A43*TY_B25*TY_B32 - TY_A43*TY_B23*TY_B34 - TY_A42*TY_B24*TY_B34 - TY_A41*TY_B25*TY_B34 + TY_A52*pow(TY_B24,2)
910        TY_F39 = 2*TY_A52*TY_B24*TY_B25 - TY_A43*TY_B24*TY_B34 - TY_A42*TY_B25*TY_B34
911        TY_F310 = -(TY_A43*TY_B25*TY_B34) + TY_A52*pow(TY_B25,2)
912       
913        if( prnt )
914                printf "\rF = \r"
915//              printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0., 0.,  0.,  TY_F14, 0.,  TY_F16, 0.,  TY_F18, 0.,  0.
916//              printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0., 0.,  TY_F23, TY_F24, TY_F25, TY_F26, TY_F27, TY_F28, TY_F29, 0. 
917//              printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0., TY_F32, TY_F33, TY_F34, TY_F35, TY_F36, TY_F37, TY_F38, TY_F39, TY_F310
918                printf "TY_F14 = %15.12g\r",TY_F14
919                printf "TY_F16 = %15.12g\r",TY_F16
920                printf "TY_F18 = %15.12g\r",TY_F18
921                printf "TY_F23 = %15.12g\r",TY_F23
922                printf "TY_F24 = %15.12g\r",TY_F24
923                printf "TY_F25 = %15.12g\r",TY_F25
924                printf "TY_F26 = %15.12g\r",TY_F26
925                printf "TY_F27 = %15.12g\r",TY_F27
926                printf "TY_F28 = %15.12g\r",TY_F28
927                printf "TY_F29 = %15.12g\r",TY_F29
928                printf "TY_F32 = %15.12g\r",TY_F32
929                printf "TY_F33 = %15.12g\r",TY_F33
930                printf "TY_F34 = %15.12g\r",TY_F34
931                printf "TY_F35 = %15.12g\r",TY_F35
932                printf "TY_F36 = %15.12g\r",TY_F36
933                printf "TY_F37 = %15.12g\r",TY_F37
934                printf "TY_F38 = %15.12g\r",TY_F38
935                printf "TY_F39 = %15.12g\r",TY_F39
936                printf "TY_F310 = %15.12g\r",TY_F310
937        endif
938       
939        NVAR TY_G13  = root:yuk:TY_G13
940        NVAR TY_G14  = root:yuk:TY_G14
941        NVAR TY_G15  = root:yuk:TY_G15
942        NVAR TY_G16  = root:yuk:TY_G16
943        NVAR TY_G17  = root:yuk:TY_G17
944        NVAR TY_G18  = root:yuk:TY_G18
945        NVAR TY_G19  = root:yuk:TY_G19
946        NVAR TY_G110 = root:yuk:TY_G110
947        NVAR TY_G111 = root:yuk:TY_G111
948        NVAR TY_G112 = root:yuk:TY_G112
949        NVAR TY_G113 = root:yuk:TY_G113
950        NVAR TY_G22  = root:yuk:TY_G22
951        NVAR TY_G23  = root:yuk:TY_G23
952        NVAR TY_G24  = root:yuk:TY_G24
953        NVAR TY_G25  = root:yuk:TY_G25
954        NVAR TY_G26  = root:yuk:TY_G26
955        NVAR TY_G27  = root:yuk:TY_G27
956        NVAR TY_G28  = root:yuk:TY_G28
957        NVAR TY_G29  = root:yuk:TY_G29
958        NVAR TY_G210 = root:yuk:TY_G210
959        NVAR TY_G211 = root:yuk:TY_G211
960        NVAR TY_G212 = root:yuk:TY_G212
961        NVAR TY_G213 = root:yuk:TY_G213
962        NVAR TY_G214 = root:yuk:TY_G214
963       
964       
965        TY_G13  = -(TY_B12*TY_F32)
966        TY_G14  = -(TY_B12*TY_F33)
967        TY_G15  = TY_B32*TY_F14 - TY_B14*TY_F32 - TY_B12*TY_F34
968        TY_G16  = -(TY_B14*TY_F33) - TY_B12*TY_F35
969        TY_G17  = TY_B34*TY_F14 + TY_B32*TY_F16 - TY_B14*TY_F34 - TY_B12*TY_F36
970        TY_G18  = -(TY_B14*TY_F35) - TY_B12*TY_F37
971        TY_G19  = TY_B34*TY_F16 + TY_B32*TY_F18 - TY_B14*TY_F36 - TY_B12*TY_F38
972        TY_G110 = -(TY_B14*TY_F37) - TY_B12*TY_F39
973        TY_G111 = TY_B34*TY_F18 - TY_B12*TY_F310 - TY_B14*TY_F38
974        TY_G112 = -(TY_B14*TY_F39)
975        TY_G113 = -(TY_B14*TY_F310)
976        TY_G22  = -(TY_B21*TY_F32)
977        TY_G23  = -(TY_B22*TY_F32) - TY_B21*TY_F33
978        TY_G24  = TY_B32*TY_F23 - TY_B23*TY_F32 - TY_B22*TY_F33 - TY_B21*TY_F34
979        TY_G25  = TY_B32*TY_F24 - TY_B24*TY_F32 - TY_B23*TY_F33 - TY_B22*TY_F34 - TY_B21*TY_F35
980        TY_G26  = TY_B34*TY_F23 + TY_B32*TY_F25 - TY_B25*TY_F32 - TY_B24*TY_F33 - TY_B23*TY_F34 - TY_B22*TY_F35 - TY_B21*TY_F36
981        TY_G27  = TY_B34*TY_F24 + TY_B32*TY_F26 - TY_B25*TY_F33 - TY_B24*TY_F34 - TY_B23*TY_F35 - TY_B22*TY_F36 - TY_B21*TY_F37
982        TY_G28  = TY_B34*TY_F25 + TY_B32*TY_F27 - TY_B25*TY_F34 - TY_B24*TY_F35 - TY_B23*TY_F36 - TY_B22*TY_F37 - TY_B21*TY_F38
983        TY_G29  = TY_B34*TY_F26 + TY_B32*TY_F28 - TY_B25*TY_F35 - TY_B24*TY_F36 - TY_B23*TY_F37 - TY_B22*TY_F38 - TY_B21*TY_F39
984        TY_G210 = TY_B34*TY_F27 + TY_B32*TY_F29 - TY_B21*TY_F310 - TY_B25*TY_F36 - TY_B24*TY_F37 - TY_B23*TY_F38 - TY_B22*TY_F39
985        TY_G211 = TY_B34*TY_F28 - TY_B22*TY_F310 - TY_B25*TY_F37 - TY_B24*TY_F38 - TY_B23*TY_F39
986        TY_G212 = TY_B34*TY_F29 - TY_B23*TY_F310 - TY_B25*TY_F38 - TY_B24*TY_F39
987        TY_G213 = -(TY_B24*TY_F310) - TY_B25*TY_F39
988        TY_G214 = -(TY_B25*TY_F310)
989       
990        if( prnt )
991                printf "\rG = \r"
992//              printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0., 0.,  TY_G13, TY_G14, TY_G15, TY_G16, TY_G17, TY_G18, TY_G19, TY_G110, TY_G111, TY_G112, TY_G113, 0.
993//              printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0., TY_G22, TY_G23, TY_G24, TY_G25, TY_G26, TY_G27, TY_G28, TY_G29, TY_G210, TY_G211, TY_G212, TY_G213, TY_G214
994                printf "TY_G13  = %15.12g\r",TY_G13
995                printf "TY_G14  = %15.12g\r",TY_G14
996                printf "TY_G15  = %15.12g\r",TY_G15
997                printf "TY_G16  = %15.12g\r",TY_G16
998                printf "TY_G17  = %15.12g\r",TY_G17
999                printf "TY_G18  = %15.12g\r",TY_G18
1000                printf "TY_G19  = %15.12g\r",TY_G19
1001                printf "TY_G110 = %15.12g\r",TY_G110
1002                printf "TY_G111 = %15.12g\r",TY_G111
1003                printf "TY_G112 = %15.12g\r",TY_G112
1004                printf "TY_G113 = %15.12g\r",TY_G113
1005                printf "TY_G22  = %15.12g\r",TY_G22
1006                printf "TY_G23  = %15.12g\r",TY_G23
1007                printf "TY_G24  = %15.12g\r",TY_G24
1008                printf "TY_G25  = %15.12g\r",TY_G25
1009                printf "TY_G26  = %15.12g\r",TY_G26
1010                printf "TY_G27  = %15.12g\r",TY_G27
1011                printf "TY_G28  = %15.12g\r",TY_G28
1012                printf "TY_G29  = %15.12g\r",TY_G29
1013                printf "TY_G210 = %15.12g\r",TY_G210
1014                printf "TY_G211 = %15.12g\r",TY_G211
1015                printf "TY_G212 = %15.12g\r",TY_G212
1016                printf "TY_G213 = %15.12g\r",TY_G213
1017                printf "TY_G214 = %15.12g\r",TY_G214
1018        endif
1019       
1020        Make/O/D/N=23 TY_w
1021       
1022        // coefficients for polynomial
1023        TY_w[0] = (-(TY_A21*TY_B12) + TY_A12*TY_B21)*(TY_A52*TY_B21 - TY_A41*TY_B32)*pow(TY_B21,2)*pow(TY_B32,3)
1024       
1025        TY_w[1] = 2*TY_B32*TY_G13*TY_G14 - TY_B24*TY_G13*TY_G22 - TY_B23*TY_G14*TY_G22 - TY_B22*TY_G15*TY_G22 - TY_B21*TY_G16*TY_G22 - TY_B23*TY_G13*TY_G23 - TY_B22*TY_G14*TY_G23
1026        TY_w[1] += - TY_B21*TY_G15*TY_G23 + 2*TY_B14*TY_G22*TY_G23 - TY_B22*TY_G13*TY_G24 - TY_B21*TY_G14*TY_G24 + 2*TY_B12*TY_G23*TY_G24 - TY_B21*TY_G13*TY_G25 + 2*TY_B12*TY_G22*TY_G25
1027       
1028        TY_w[2] = -(TY_B25*TY_G13*TY_G22) - TY_B24*TY_G14*TY_G22 - TY_B23*TY_G15*TY_G22 - TY_B22*TY_G16*TY_G22 - TY_B21*TY_G17*TY_G22 - TY_B24*TY_G13*TY_G23 - TY_B23*TY_G14*TY_G23 - TY_B22*TY_G15*TY_G23 - TY_B21*TY_G16*TY_G23
1029        TY_w[2] += -TY_B23*TY_G13*TY_G24 - TY_B22*TY_G14*TY_G24 - TY_B21*TY_G15*TY_G24 + 2*TY_B14*TY_G22*TY_G24 - TY_B22*TY_G13*TY_G25 - TY_B21*TY_G14*TY_G25 + 2*TY_B12*TY_G23*TY_G25 - TY_B21*TY_G13*TY_G26 + 2*TY_B12*TY_G22*TY_G26
1030        TY_w[2] += +TY_B34*pow(TY_G13,2) + TY_B32*(2*TY_G13*TY_G15 + pow(TY_G14,2)) + TY_B14*pow(TY_G23,2) + TY_B12*pow(TY_G24,2)
1031       
1032        TY_w[3] = 2*TY_B34*TY_G13*TY_G14 + 2*TY_B32*(TY_G14*TY_G15 + TY_G13*TY_G16) - TY_B25*TY_G14*TY_G22 - TY_B24*TY_G15*TY_G22 - TY_B23*TY_G16*TY_G22 - TY_B22*TY_G17*TY_G22 - TY_B21*TY_G18*TY_G22 - TY_B25*TY_G13*TY_G23 
1033        TY_w[3] += -TY_B24*TY_G14*TY_G23 - TY_B23*TY_G15*TY_G23 - TY_B22*TY_G16*TY_G23 - TY_B21*TY_G17*TY_G23 - TY_B24*TY_G13*TY_G24 - TY_B23*TY_G14*TY_G24 - TY_B22*TY_G15*TY_G24 - TY_B21*TY_G16*TY_G24 + 2*TY_B14*TY_G23*TY_G24 
1034        TY_w[3] += -TY_B23*TY_G13*TY_G25 - TY_B22*TY_G14*TY_G25 - TY_B21*TY_G15*TY_G25 + 2*TY_B14*TY_G22*TY_G25 + 2*TY_B12*TY_G24*TY_G25 - TY_B22*TY_G13*TY_G26 - TY_B21*TY_G14*TY_G26 + 2*TY_B12*TY_G23*TY_G26 - TY_B21*TY_G13*TY_G27
1035        TY_w[3] += 2*TY_B12*TY_G22*TY_G27;
1036       
1037        TY_w[4] = -(TY_B25*TY_G15*TY_G22) - TY_B24*TY_G16*TY_G22 - TY_B23*TY_G17*TY_G22 - TY_B22*TY_G18*TY_G22 - TY_B21*TY_G19*TY_G22 - TY_B25*TY_G14*TY_G23 - TY_B24*TY_G15*TY_G23 - TY_B23*TY_G16*TY_G23 - TY_B22*TY_G17*TY_G23 
1038        TY_w[4] += -TY_B21*TY_G18*TY_G23 - TY_B25*TY_G13*TY_G24 - TY_B24*TY_G14*TY_G24 - TY_B23*TY_G15*TY_G24 - TY_B22*TY_G16*TY_G24 - TY_B21*TY_G17*TY_G24 - TY_B24*TY_G13*TY_G25 - TY_B23*TY_G14*TY_G25 - TY_B22*TY_G15*TY_G25 
1039        TY_w[4] += -TY_B21*TY_G16*TY_G25 + 2*TY_B14*TY_G23*TY_G25 - TY_B23*TY_G13*TY_G26 - TY_B22*TY_G14*TY_G26 - TY_B21*TY_G15*TY_G26 + 2*TY_B14*TY_G22*TY_G26 + 2*TY_B12*TY_G24*TY_G26 - TY_B22*TY_G13*TY_G27 - TY_B21*TY_G14*TY_G27 
1040        TY_w[4] += 2*TY_B12*TY_G23*TY_G27 - TY_B21*TY_G13*TY_G28 + 2*TY_B12*TY_G22*TY_G28 + TY_B34*(2*TY_G13*TY_G15 + pow(TY_G14,2)) + TY_B32*(2*TY_G14*TY_G16 + 2*TY_G13*TY_G17 + pow(TY_G15,2)) + TY_B14*pow(TY_G24,2) 
1041        TY_w[4] += TY_B12*pow(TY_G25,2)
1042       
1043        TY_w[5] = 2*TY_B34*(TY_G14*TY_G15 + TY_G13*TY_G16) + 2*TY_B32*(TY_G15*TY_G16 + TY_G14*TY_G17 + TY_G13*TY_G18) - TY_B21*TY_G110*TY_G22 - TY_B25*TY_G16*TY_G22 - TY_B24*TY_G17*TY_G22 - TY_B23*TY_G18*TY_G22 - TY_B22*TY_G19*TY_G22
1044        TY_w[5] += -TY_B25*TY_G15*TY_G23 - TY_B24*TY_G16*TY_G23 - TY_B23*TY_G17*TY_G23 - TY_B22*TY_G18*TY_G23 - TY_B21*TY_G19*TY_G23 - TY_B25*TY_G14*TY_G24 - TY_B24*TY_G15*TY_G24 - TY_B23*TY_G16*TY_G24 - TY_B22*TY_G17*TY_G24
1045        TY_w[5] += -TY_B21*TY_G18*TY_G24 - TY_B25*TY_G13*TY_G25 - TY_B24*TY_G14*TY_G25 - TY_B23*TY_G15*TY_G25 - TY_B22*TY_G16*TY_G25 - TY_B21*TY_G17*TY_G25 + 2*TY_B14*TY_G24*TY_G25 - TY_B24*TY_G13*TY_G26 - TY_B23*TY_G14*TY_G26 
1046        TY_w[5] += -TY_B22*TY_G15*TY_G26 - TY_B21*TY_G16*TY_G26 + 2*TY_B14*TY_G23*TY_G26 + 2*TY_B12*TY_G25*TY_G26 - TY_B23*TY_G13*TY_G27 - TY_B22*TY_G14*TY_G27 - TY_B21*TY_G15*TY_G27 + 2*TY_B14*TY_G22*TY_G27 + 2*TY_B12*TY_G24*TY_G27 
1047        TY_w[5] += -TY_B22*TY_G13*TY_G28 - TY_B21*TY_G14*TY_G28 + 2*TY_B12*TY_G23*TY_G28 - TY_B21*TY_G13*TY_G29 + 2*TY_B12*TY_G22*TY_G29
1048       
1049        TY_w[6] = -(TY_B22*TY_G110*TY_G22) - TY_B21*TY_G111*TY_G22 - TY_B25*TY_G17*TY_G22 - TY_B24*TY_G18*TY_G22 - TY_B23*TY_G19*TY_G22 + TY_G210*(-(TY_B21*TY_G13) + 2*TY_B12*TY_G22) - TY_B21*TY_G110*TY_G23 - TY_B25*TY_G16*TY_G23 
1050        TY_w[6] += -TY_B24*TY_G17*TY_G23 - TY_B23*TY_G18*TY_G23 - TY_B22*TY_G19*TY_G23 - TY_B25*TY_G15*TY_G24 - TY_B24*TY_G16*TY_G24 - TY_B23*TY_G17*TY_G24 - TY_B22*TY_G18*TY_G24 - TY_B21*TY_G19*TY_G24 - TY_B25*TY_G14*TY_G25 
1051        TY_w[6] += -TY_B24*TY_G15*TY_G25 - TY_B23*TY_G16*TY_G25 - TY_B22*TY_G17*TY_G25 - TY_B21*TY_G18*TY_G25 - TY_B25*TY_G13*TY_G26 - TY_B24*TY_G14*TY_G26 - TY_B23*TY_G15*TY_G26 - TY_B22*TY_G16*TY_G26 - TY_B21*TY_G17*TY_G26 
1052        TY_w[6] += 2*TY_B14*TY_G24*TY_G26 - TY_B24*TY_G13*TY_G27 - TY_B23*TY_G14*TY_G27 - TY_B22*TY_G15*TY_G27 - TY_B21*TY_G16*TY_G27 + 2*TY_B14*TY_G23*TY_G27 + 2*TY_B12*TY_G25*TY_G27 - TY_B23*TY_G13*TY_G28 - TY_B22*TY_G14*TY_G28 
1053        TY_w[6] += -TY_B21*TY_G15*TY_G28 + 2*TY_B14*TY_G22*TY_G28 + 2*TY_B12*TY_G24*TY_G28 - TY_B22*TY_G13*TY_G29 - TY_B21*TY_G14*TY_G29 + 2*TY_B12*TY_G23*TY_G29 + TY_B34*(2*TY_G14*TY_G16 + 2*TY_G13*TY_G17 + pow(TY_G15,2)) 
1054        TY_w[6] += TY_B32*(2*(TY_G15*TY_G17 + TY_G14*TY_G18 + TY_G13*TY_G19) + pow(TY_G16,2)) + TY_B14*pow(TY_G25,2) + TY_B12*pow(TY_G26,2)
1055       
1056        TY_w[7] = 2*TY_B34*(TY_G15*TY_G16 + TY_G14*TY_G17 + TY_G13*TY_G18) + 2*TY_B32*(TY_G110*TY_G13 + TY_G16*TY_G17 + TY_G15*TY_G18 + TY_G14*TY_G19) - TY_B22*TY_G13*TY_G210 - TY_B21*TY_G14*TY_G210 - TY_B23*TY_G110*TY_G22 
1057        TY_w[7] += -TY_B22*TY_G111*TY_G22 - TY_B21*TY_G112*TY_G22 - TY_B25*TY_G18*TY_G22 - TY_B24*TY_G19*TY_G22 + TY_G211*(-(TY_B21*TY_G13) + 2*TY_B12*TY_G22) - TY_B22*TY_G110*TY_G23 - TY_B21*TY_G111*TY_G23 - TY_B25*TY_G17*TY_G23 
1058        TY_w[7] += -TY_B24*TY_G18*TY_G23 - TY_B23*TY_G19*TY_G23 + 2*TY_B12*TY_G210*TY_G23 - TY_B21*TY_G110*TY_G24 - TY_B25*TY_G16*TY_G24 - TY_B24*TY_G17*TY_G24 - TY_B23*TY_G18*TY_G24 - TY_B22*TY_G19*TY_G24 - TY_B25*TY_G15*TY_G25 
1059        TY_w[7] += -TY_B24*TY_G16*TY_G25 - TY_B23*TY_G17*TY_G25 - TY_B22*TY_G18*TY_G25 - TY_B21*TY_G19*TY_G25 - TY_B25*TY_G14*TY_G26 - TY_B24*TY_G15*TY_G26 - TY_B23*TY_G16*TY_G26 - TY_B22*TY_G17*TY_G26 - TY_B21*TY_G18*TY_G26
1060        TY_w[7] += 2*TY_B14*TY_G25*TY_G26 - TY_B25*TY_G13*TY_G27 - TY_B24*TY_G14*TY_G27 - TY_B23*TY_G15*TY_G27 - TY_B22*TY_G16*TY_G27 - TY_B21*TY_G17*TY_G27 + 2*TY_B14*TY_G24*TY_G27 + 2*TY_B12*TY_G26*TY_G27 - TY_B24*TY_G13*TY_G28 
1061        TY_w[7] += -TY_B23*TY_G14*TY_G28 - TY_B22*TY_G15*TY_G28 - TY_B21*TY_G16*TY_G28 + 2*TY_B14*TY_G23*TY_G28 + 2*TY_B12*TY_G25*TY_G28 - TY_B23*TY_G13*TY_G29 - TY_B22*TY_G14*TY_G29 - TY_B21*TY_G15*TY_G29 + 2*TY_B14*TY_G22*TY_G29 
1062        TY_w[7] += 2*TY_B12*TY_G24*TY_G29
1063       
1064        TY_w[8] = -(TY_B23*TY_G13*TY_G210) - TY_B22*TY_G14*TY_G210 - TY_B21*TY_G15*TY_G210 - TY_B22*TY_G13*TY_G211 - TY_B21*TY_G14*TY_G211 - TY_B21*TY_G13*TY_G212 - TY_B24*TY_G110*TY_G22 - TY_B23*TY_G111*TY_G22 - TY_B22*TY_G112*TY_G22 
1065        TY_w[8] += -TY_B21*TY_G113*TY_G22 - TY_B25*TY_G19*TY_G22 + 2*TY_B14*TY_G210*TY_G22 + 2*TY_B12*TY_G212*TY_G22 - TY_B23*TY_G110*TY_G23 - TY_B22*TY_G111*TY_G23 - TY_B21*TY_G112*TY_G23 - TY_B25*TY_G18*TY_G23 - TY_B24*TY_G19*TY_G23 
1066        TY_w[8] += 2*TY_B12*TY_G211*TY_G23 - TY_B22*TY_G110*TY_G24 - TY_B21*TY_G111*TY_G24 - TY_B25*TY_G17*TY_G24 - TY_B24*TY_G18*TY_G24 - TY_B23*TY_G19*TY_G24 + 2*TY_B12*TY_G210*TY_G24 - TY_B21*TY_G110*TY_G25 - TY_B25*TY_G16*TY_G25 
1067        TY_w[8] += -TY_B24*TY_G17*TY_G25 - TY_B23*TY_G18*TY_G25 - TY_B22*TY_G19*TY_G25 - TY_B25*TY_G15*TY_G26 - TY_B24*TY_G16*TY_G26 - TY_B23*TY_G17*TY_G26 - TY_B22*TY_G18*TY_G26 - TY_B21*TY_G19*TY_G26 - TY_B25*TY_G14*TY_G27
1068        TY_w[8] += -TY_B24*TY_G15*TY_G27 - TY_B23*TY_G16*TY_G27 - TY_B22*TY_G17*TY_G27 - TY_B21*TY_G18*TY_G27 + 2*TY_B14*TY_G25*TY_G27 - TY_B25*TY_G13*TY_G28 - TY_B24*TY_G14*TY_G28 - TY_B23*TY_G15*TY_G28 - TY_B22*TY_G16*TY_G28 
1069        TY_w[8] += -TY_B21*TY_G17*TY_G28 + 2*TY_B14*TY_G24*TY_G28 + 2*TY_B12*TY_G26*TY_G28 - TY_B24*TY_G13*TY_G29 - TY_B23*TY_G14*TY_G29 - TY_B22*TY_G15*TY_G29 - TY_B21*TY_G16*TY_G29 + 2*TY_B14*TY_G23*TY_G29 + 2*TY_B12*TY_G25*TY_G29 
1070        TY_w[8] += TY_B34*(2*(TY_G15*TY_G17 + TY_G14*TY_G18 + TY_G13*TY_G19) + pow(TY_G16,2)) + TY_B32*(2*(TY_G111*TY_G13 + TY_G110*TY_G14 + TY_G16*TY_G18 + TY_G15*TY_G19) + pow(TY_G17,2)) + TY_B14*pow(TY_G26,2)
1071        TY_w[8] += TY_B12*pow(TY_G27,2)
1072       
1073        TY_w[9] = 2*TY_B34*(TY_G110*TY_G13 + TY_G16*TY_G17 + TY_G15*TY_G18 + TY_G14*TY_G19) + 2*TY_B32*(TY_G112*TY_G13 + TY_G111*TY_G14 + TY_G110*TY_G15 + TY_G17*TY_G18 + TY_G16*TY_G19) - TY_B24*TY_G13*TY_G210 - TY_B23*TY_G14*TY_G210 
1074        TY_w[9] += -TY_B22*TY_G15*TY_G210 - TY_B21*TY_G16*TY_G210 - TY_B23*TY_G13*TY_G211 - TY_B22*TY_G14*TY_G211 - TY_B21*TY_G15*TY_G211 - TY_B22*TY_G13*TY_G212 - TY_B21*TY_G14*TY_G212 - TY_B25*TY_G110*TY_G22 - TY_B24*TY_G111*TY_G22 
1075        TY_w[9] += -TY_B23*TY_G112*TY_G22 - TY_B22*TY_G113*TY_G22 + 2*TY_B14*TY_G211*TY_G22 + TY_G213*(-(TY_B21*TY_G13) + 2*TY_B12*TY_G22) - TY_B24*TY_G110*TY_G23 - TY_B23*TY_G111*TY_G23 - TY_B22*TY_G112*TY_G23
1076        TY_w[9] += -TY_B21*TY_G113*TY_G23 - TY_B25*TY_G19*TY_G23 + 2*TY_B14*TY_G210*TY_G23 + 2*TY_B12*TY_G212*TY_G23 - TY_B23*TY_G110*TY_G24 - TY_B22*TY_G111*TY_G24 - TY_B21*TY_G112*TY_G24 - TY_B25*TY_G18*TY_G24 - TY_B24*TY_G19*TY_G24 
1077        TY_w[9] += 2*TY_B12*TY_G211*TY_G24 - TY_B22*TY_G110*TY_G25 - TY_B21*TY_G111*TY_G25 - TY_B25*TY_G17*TY_G25 - TY_B24*TY_G18*TY_G25 - TY_B23*TY_G19*TY_G25 + 2*TY_B12*TY_G210*TY_G25 - TY_B21*TY_G110*TY_G26 - TY_B25*TY_G16*TY_G26 
1078        TY_w[9] += -TY_B24*TY_G17*TY_G26 - TY_B23*TY_G18*TY_G26 - TY_B22*TY_G19*TY_G26 - TY_B25*TY_G15*TY_G27 - TY_B24*TY_G16*TY_G27 - TY_B23*TY_G17*TY_G27 - TY_B22*TY_G18*TY_G27 - TY_B21*TY_G19*TY_G27 + 2*TY_B14*TY_G26*TY_G27
1079        TY_w[9] += -TY_B25*TY_G14*TY_G28 - TY_B24*TY_G15*TY_G28 - TY_B23*TY_G16*TY_G28 - TY_B22*TY_G17*TY_G28 - TY_B21*TY_G18*TY_G28 + 2*TY_B14*TY_G25*TY_G28 + 2*TY_B12*TY_G27*TY_G28 - TY_B25*TY_G13*TY_G29 - TY_B24*TY_G14*TY_G29 
1080        TY_w[9] += -TY_B23*TY_G15*TY_G29 - TY_B22*TY_G16*TY_G29 - TY_B21*TY_G17*TY_G29 + 2*TY_B14*TY_G24*TY_G29 + 2*TY_B12*TY_G26*TY_G29
1081       
1082        TY_w[10] = -(TY_B25*TY_G13*TY_G210) - TY_B24*TY_G14*TY_G210 - TY_B23*TY_G15*TY_G210 - TY_B22*TY_G16*TY_G210 - TY_B21*TY_G17*TY_G210 - TY_B24*TY_G13*TY_G211 - TY_B23*TY_G14*TY_G211 - TY_B22*TY_G15*TY_G211 - TY_B21*TY_G16*TY_G211 
1083        TY_w[10] += -TY_B23*TY_G13*TY_G212 - TY_B22*TY_G14*TY_G212 - TY_B21*TY_G15*TY_G212 - TY_B22*TY_G13*TY_G213 - TY_B21*TY_G14*TY_G213 - TY_B21*TY_G13*TY_G214 - TY_B25*TY_G111*TY_G22 - TY_B24*TY_G112*TY_G22 - TY_B23*TY_G113*TY_G22 
1084        TY_w[10] += 2*TY_B14*TY_G212*TY_G22 + 2*TY_B12*TY_G214*TY_G22 - TY_B25*TY_G110*TY_G23 - TY_B24*TY_G111*TY_G23 - TY_B23*TY_G112*TY_G23 - TY_B22*TY_G113*TY_G23 + 2*TY_B14*TY_G211*TY_G23 + 2*TY_B12*TY_G213*TY_G23
1085        TY_w[10] += -TY_B24*TY_G110*TY_G24 - TY_B23*TY_G111*TY_G24 - TY_B22*TY_G112*TY_G24 - TY_B21*TY_G113*TY_G24 - TY_B25*TY_G19*TY_G24 + 2*TY_B14*TY_G210*TY_G24 + 2*TY_B12*TY_G212*TY_G24 - TY_B23*TY_G110*TY_G25
1086        TY_w[10] += -TY_B22*TY_G111*TY_G25 - TY_B21*TY_G112*TY_G25 - TY_B25*TY_G18*TY_G25 - TY_B24*TY_G19*TY_G25 + 2*TY_B12*TY_G211*TY_G25 - TY_B22*TY_G110*TY_G26 - TY_B21*TY_G111*TY_G26 - TY_B25*TY_G17*TY_G26 - TY_B24*TY_G18*TY_G26 
1087        TY_w[10] += -TY_B23*TY_G19*TY_G26 + 2*TY_B12*TY_G210*TY_G26 - TY_B21*TY_G110*TY_G27 - TY_B25*TY_G16*TY_G27 - TY_B24*TY_G17*TY_G27 - TY_B23*TY_G18*TY_G27 - TY_B22*TY_G19*TY_G27 - TY_B25*TY_G15*TY_G28 - TY_B24*TY_G16*TY_G28
1088        TY_w[10] += -TY_B23*TY_G17*TY_G28 - TY_B22*TY_G18*TY_G28 - TY_B21*TY_G19*TY_G28 + 2*TY_B14*TY_G26*TY_G28 - TY_B25*TY_G14*TY_G29 - TY_B24*TY_G15*TY_G29 - TY_B23*TY_G16*TY_G29 - TY_B22*TY_G17*TY_G29 - TY_B21*TY_G18*TY_G29
1089        TY_w[10] += 2*TY_B14*TY_G25*TY_G29 + 2*TY_B12*TY_G27*TY_G29 + TY_B34*(2*(TY_G111*TY_G13 + TY_G110*TY_G14 + TY_G16*TY_G18 + TY_G15*TY_G19) + pow(TY_G17,2))
1090        TY_w[10] += TY_B32*(2*(TY_G113*TY_G13 + TY_G112*TY_G14 + TY_G111*TY_G15 + TY_G110*TY_G16 + TY_G17*TY_G19) + pow(TY_G18,2)) + TY_B14*pow(TY_G27,2) + TY_B12*pow(TY_G28,2)
1091       
1092        TY_w[11] = 2*TY_B34*(TY_G112*TY_G13 + TY_G111*TY_G14 + TY_G110*TY_G15 + TY_G17*TY_G18 + TY_G16*TY_G19) + 2*TY_B32*(TY_G113*TY_G14 + TY_G112*TY_G15 + TY_G111*TY_G16 + TY_G110*TY_G17 + TY_G18*TY_G19) - TY_B25*TY_G14*TY_G210 
1093        TY_w[11] += -TY_B24*TY_G15*TY_G210 - TY_B23*TY_G16*TY_G210 - TY_B22*TY_G17*TY_G210 - TY_B21*TY_G18*TY_G210 - TY_B25*TY_G13*TY_G211 - TY_B24*TY_G14*TY_G211 - TY_B23*TY_G15*TY_G211 - TY_B22*TY_G16*TY_G211 - TY_B21*TY_G17*TY_G211 
1094        TY_w[11] += -TY_B24*TY_G13*TY_G212 - TY_B23*TY_G14*TY_G212 - TY_B22*TY_G15*TY_G212 - TY_B21*TY_G16*TY_G212 - TY_B23*TY_G13*TY_G213 - TY_B22*TY_G14*TY_G213 - TY_B21*TY_G15*TY_G213 - TY_B25*TY_G112*TY_G22 - TY_B24*TY_G113*TY_G22 
1095        TY_w[11] += 2*TY_B14*TY_G213*TY_G22 - TY_B25*TY_G111*TY_G23 - TY_B24*TY_G112*TY_G23 - TY_B23*TY_G113*TY_G23 + 2*TY_B14*TY_G212*TY_G23 - TY_G214*(TY_B22*TY_G13 + TY_B21*TY_G14 - 2*TY_B12*TY_G23) - TY_B25*TY_G110*TY_G24
1096        TY_w[11] += -TY_B24*TY_G111*TY_G24 - TY_B23*TY_G112*TY_G24 - TY_B22*TY_G113*TY_G24 + 2*TY_B14*TY_G211*TY_G24 + 2*TY_B12*TY_G213*TY_G24 - TY_B24*TY_G110*TY_G25 - TY_B23*TY_G111*TY_G25 - TY_B22*TY_G112*TY_G25
1097        TY_w[11] += -TY_B21*TY_G113*TY_G25 - TY_B25*TY_G19*TY_G25 + 2*TY_B14*TY_G210*TY_G25 + 2*TY_B12*TY_G212*TY_G25 - TY_B23*TY_G110*TY_G26 - TY_B22*TY_G111*TY_G26 - TY_B21*TY_G112*TY_G26 - TY_B25*TY_G18*TY_G26 - TY_B24*TY_G19*TY_G26 
1098        TY_w[11] += 2*TY_B12*TY_G211*TY_G26 - TY_B22*TY_G110*TY_G27 - TY_B21*TY_G111*TY_G27 - TY_B25*TY_G17*TY_G27 - TY_B24*TY_G18*TY_G27 - TY_B23*TY_G19*TY_G27 + 2*TY_B12*TY_G210*TY_G27 - TY_B21*TY_G110*TY_G28 - TY_B25*TY_G16*TY_G28
1099        TY_w[11] += -TY_B24*TY_G17*TY_G28 - TY_B23*TY_G18*TY_G28 - TY_B22*TY_G19*TY_G28 + 2*TY_B14*TY_G27*TY_G28 - TY_B25*TY_G15*TY_G29 - TY_B24*TY_G16*TY_G29 - TY_B23*TY_G17*TY_G29 - TY_B22*TY_G18*TY_G29 - TY_B21*TY_G19*TY_G29
1100        TY_w[11] += 2*TY_B14*TY_G26*TY_G29 + 2*TY_B12*TY_G28*TY_G29
1101       
1102        TY_w[12] = -(TY_B25*TY_G15*TY_G210) - TY_B24*TY_G16*TY_G210 - TY_B23*TY_G17*TY_G210 - TY_B22*TY_G18*TY_G210 - TY_B21*TY_G19*TY_G210 - TY_B25*TY_G14*TY_G211 - TY_B24*TY_G15*TY_G211 - TY_B23*TY_G16*TY_G211 - TY_B22*TY_G17*TY_G211
1103        TY_w[12] += -TY_B21*TY_G18*TY_G211 - TY_B25*TY_G13*TY_G212 - TY_B24*TY_G14*TY_G212 - TY_B23*TY_G15*TY_G212 - TY_B22*TY_G16*TY_G212 - TY_B21*TY_G17*TY_G212 - TY_B24*TY_G13*TY_G213 - TY_B23*TY_G14*TY_G213 - TY_B22*TY_G15*TY_G213 
1104        TY_w[12] += -TY_B21*TY_G16*TY_G213 - TY_B25*TY_G113*TY_G22 - TY_B25*TY_G112*TY_G23 - TY_B24*TY_G113*TY_G23 + 2*TY_B14*TY_G213*TY_G23 - TY_B25*TY_G111*TY_G24 - TY_B24*TY_G112*TY_G24 - TY_B23*TY_G113*TY_G24
1105        TY_w[12] += 2*TY_B14*TY_G212*TY_G24 - TY_G214*(TY_B23*TY_G13 + TY_B22*TY_G14 + TY_B21*TY_G15 - 2*TY_B14*TY_G22 - 2*TY_B12*TY_G24) - TY_B25*TY_G110*TY_G25 - TY_B24*TY_G111*TY_G25 - TY_B23*TY_G112*TY_G25
1106        TY_w[12] += -TY_B22*TY_G113*TY_G25 + 2*TY_B14*TY_G211*TY_G25 + 2*TY_B12*TY_G213*TY_G25 - TY_B24*TY_G110*TY_G26 - TY_B23*TY_G111*TY_G26 - TY_B22*TY_G112*TY_G26 - TY_B21*TY_G113*TY_G26 - TY_B25*TY_G19*TY_G26 
1107        TY_w[12] += 2*TY_B14*TY_G210*TY_G26 + 2*TY_B12*TY_G212*TY_G26 - TY_B23*TY_G110*TY_G27 - TY_B22*TY_G111*TY_G27 - TY_B21*TY_G112*TY_G27 - TY_B25*TY_G18*TY_G27 - TY_B24*TY_G19*TY_G27 + 2*TY_B12*TY_G211*TY_G27 
1108        TY_w[12] += -TY_B22*TY_G110*TY_G28 - TY_B21*TY_G111*TY_G28 - TY_B25*TY_G17*TY_G28 - TY_B24*TY_G18*TY_G28 - TY_B23*TY_G19*TY_G28 + 2*TY_B12*TY_G210*TY_G28 - TY_B21*TY_G110*TY_G29 - TY_B25*TY_G16*TY_G29 - TY_B24*TY_G17*TY_G29 
1109        TY_w[12] += -TY_B23*TY_G18*TY_G29 - TY_B22*TY_G19*TY_G29 + 2*TY_B14*TY_G27*TY_G29 + TY_B34*(2*(TY_G113*TY_G13 + TY_G112*TY_G14 + TY_G111*TY_G15 + TY_G110*TY_G16 + TY_G17*TY_G19) + pow(TY_G18,2))
1110        TY_w[12] += TY_B32*(2*(TY_G113*TY_G15 + TY_G112*TY_G16 + TY_G111*TY_G17 + TY_G110*TY_G18) + pow(TY_G19,2)) + TY_B14*pow(TY_G28,2) + TY_B12*pow(TY_G29,2)
1111       
1112        TY_w[13] = 2*TY_B32*(TY_G113*TY_G16 + TY_G112*TY_G17 + TY_G111*TY_G18 + TY_G110*TY_G19) + 2*TY_B34*(TY_G113*TY_G14 + TY_G112*TY_G15 + TY_G111*TY_G16 + TY_G110*TY_G17 + TY_G18*TY_G19) - TY_B21*TY_G110*TY_G210
1113        TY_w[13] += -TY_B25*TY_G16*TY_G210 - TY_B24*TY_G17*TY_G210 - TY_B23*TY_G18*TY_G210 - TY_B22*TY_G19*TY_G210 - TY_B25*TY_G15*TY_G211 - TY_B24*TY_G16*TY_G211 - TY_B23*TY_G17*TY_G211 - TY_B22*TY_G18*TY_G211 - TY_B21*TY_G19*TY_G211 
1114        TY_w[13] += -TY_B25*TY_G14*TY_G212 - TY_B24*TY_G15*TY_G212 - TY_B23*TY_G16*TY_G212 - TY_B22*TY_G17*TY_G212 - TY_B21*TY_G18*TY_G212 - TY_B25*TY_G13*TY_G213 - TY_B24*TY_G14*TY_G213 - TY_B23*TY_G15*TY_G213 - TY_B22*TY_G16*TY_G213 
1115        TY_w[13] += -TY_B21*TY_G17*TY_G213 - TY_B25*TY_G113*TY_G23 - TY_B25*TY_G112*TY_G24 - TY_B24*TY_G113*TY_G24 + 2*TY_B14*TY_G213*TY_G24 - TY_B25*TY_G111*TY_G25 - TY_B24*TY_G112*TY_G25 - TY_B23*TY_G113*TY_G25
1116        TY_w[13] += 2*TY_B14*TY_G212*TY_G25 - TY_G214*(TY_B24*TY_G13 + TY_B23*TY_G14 + TY_B22*TY_G15 + TY_B21*TY_G16 - 2*TY_B14*TY_G23 - 2*TY_B12*TY_G25) - TY_B25*TY_G110*TY_G26 - TY_B24*TY_G111*TY_G26 - TY_B23*TY_G112*TY_G26 
1117        TY_w[13] += -TY_B22*TY_G113*TY_G26 + 2*TY_B14*TY_G211*TY_G26 + 2*TY_B12*TY_G213*TY_G26 - TY_B24*TY_G110*TY_G27 - TY_B23*TY_G111*TY_G27 - TY_B22*TY_G112*TY_G27 - TY_B21*TY_G113*TY_G27 - TY_B25*TY_G19*TY_G27
1118        TY_w[13] += 2*TY_B14*TY_G210*TY_G27 + 2*TY_B12*TY_G212*TY_G27 - TY_B23*TY_G110*TY_G28 - TY_B22*TY_G111*TY_G28 - TY_B21*TY_G112*TY_G28 - TY_B25*TY_G18*TY_G28 - TY_B24*TY_G19*TY_G28 + 2*TY_B12*TY_G211*TY_G28 
1119        TY_w[13] += -TY_B22*TY_G110*TY_G29 - TY_B21*TY_G111*TY_G29 - TY_B25*TY_G17*TY_G29 - TY_B24*TY_G18*TY_G29 - TY_B23*TY_G19*TY_G29 + 2*TY_B12*TY_G210*TY_G29 + 2*TY_B14*TY_G28*TY_G29
1120       
1121        TY_w[14] = -(TY_B22*TY_G110*TY_G210) - TY_B21*TY_G111*TY_G210 - TY_B25*TY_G17*TY_G210 - TY_B24*TY_G18*TY_G210 - TY_B23*TY_G19*TY_G210 - TY_B21*TY_G110*TY_G211 - TY_B25*TY_G16*TY_G211 - TY_B24*TY_G17*TY_G211
1122        TY_w[14] += -TY_B23*TY_G18*TY_G211 - TY_B22*TY_G19*TY_G211 - TY_B25*TY_G15*TY_G212 - TY_B24*TY_G16*TY_G212 - TY_B23*TY_G17*TY_G212 - TY_B22*TY_G18*TY_G212 - TY_B21*TY_G19*TY_G212 - TY_B25*TY_G14*TY_G213 - TY_B24*TY_G15*TY_G213
1123        TY_w[14] += -TY_B23*TY_G16*TY_G213 - TY_B22*TY_G17*TY_G213 - TY_B21*TY_G18*TY_G213 - TY_B25*TY_G113*TY_G24 - TY_B25*TY_G112*TY_G25 - TY_B24*TY_G113*TY_G25 + 2*TY_B14*TY_G213*TY_G25 - TY_B25*TY_G111*TY_G26 - TY_B24*TY_G112*TY_G26 
1124        TY_w[14] += -TY_B23*TY_G113*TY_G26 + 2*TY_B14*TY_G212*TY_G26 - TY_G214*(TY_B25*TY_G13 + TY_B24*TY_G14 + TY_B23*TY_G15 + TY_B22*TY_G16 + TY_B21*TY_G17 - 2*TY_B14*TY_G24 - 2*TY_B12*TY_G26) - TY_B25*TY_G110*TY_G27
1125        TY_w[14] += -TY_B24*TY_G111*TY_G27 - TY_B23*TY_G112*TY_G27 - TY_B22*TY_G113*TY_G27 + 2*TY_B14*TY_G211*TY_G27 + 2*TY_B12*TY_G213*TY_G27 - TY_B24*TY_G110*TY_G28 - TY_B23*TY_G111*TY_G28 - TY_B22*TY_G112*TY_G28
1126        TY_w[14] += -TY_B21*TY_G113*TY_G28 - TY_B25*TY_G19*TY_G28 + 2*TY_B14*TY_G210*TY_G28 + 2*TY_B12*TY_G212*TY_G28 - TY_B23*TY_G110*TY_G29 - TY_B22*TY_G111*TY_G29 - TY_B21*TY_G112*TY_G29 - TY_B25*TY_G18*TY_G29 - TY_B24*TY_G19*TY_G29
1127        TY_w[14] += 2*TY_B12*TY_G211*TY_G29 + TY_B32*(2*(TY_G113*TY_G17 + TY_G112*TY_G18 + TY_G111*TY_G19) + pow(TY_G110,2))
1128        TY_w[14] += TY_B34*(2*(TY_G113*TY_G15 + TY_G112*TY_G16 + TY_G111*TY_G17 + TY_G110*TY_G18) + pow(TY_G19,2)) + TY_B12*pow(TY_G210,2) + TY_B14*pow(TY_G29,2)
1129       
1130        TY_w[15] = 2*TY_B34*(TY_G113*TY_G16 + TY_G112*TY_G17 + TY_G111*TY_G18 + TY_G110*TY_G19) + 2*TY_B32*(TY_G110*TY_G111 + TY_G113*TY_G18 + TY_G112*TY_G19) - TY_B23*TY_G110*TY_G210 - TY_B22*TY_G111*TY_G210 
1131        TY_w[15] += -TY_B21*TY_G112*TY_G210 - TY_B25*TY_G18*TY_G210 - TY_B24*TY_G19*TY_G210 - TY_B22*TY_G110*TY_G211 - TY_B21*TY_G111*TY_G211 - TY_B25*TY_G17*TY_G211 - TY_B24*TY_G18*TY_G211 - TY_B23*TY_G19*TY_G211 
1132        TY_w[15] += 2*TY_B12*TY_G210*TY_G211 - TY_B21*TY_G110*TY_G212 - TY_B25*TY_G16*TY_G212 - TY_B24*TY_G17*TY_G212 - TY_B23*TY_G18*TY_G212 - TY_B22*TY_G19*TY_G212 - TY_B25*TY_G15*TY_G213 - TY_B24*TY_G16*TY_G213 
1133        TY_w[15] += -TY_B23*TY_G17*TY_G213 - TY_B22*TY_G18*TY_G213 - TY_B21*TY_G19*TY_G213 - TY_B25*TY_G113*TY_G25 - TY_B25*TY_G112*TY_G26 - TY_B24*TY_G113*TY_G26 + 2*TY_B14*TY_G213*TY_G26 - TY_B25*TY_G111*TY_G27 - TY_B24*TY_G112*TY_G27
1134        TY_w[15] += -TY_B23*TY_G113*TY_G27 + 2*TY_B14*TY_G212*TY_G27 - TY_G214*(TY_B25*TY_G14 + TY_B24*TY_G15 + TY_B23*TY_G16 + TY_B22*TY_G17 + TY_B21*TY_G18 - 2*TY_B14*TY_G25 - 2*TY_B12*TY_G27) - TY_B25*TY_G110*TY_G28
1135        TY_w[15] += -TY_B24*TY_G111*TY_G28 - TY_B23*TY_G112*TY_G28 - TY_B22*TY_G113*TY_G28 + 2*TY_B14*TY_G211*TY_G28 + 2*TY_B12*TY_G213*TY_G28 - TY_B24*TY_G110*TY_G29 - TY_B23*TY_G111*TY_G29 - TY_B22*TY_G112*TY_G29
1136        TY_w[15] += -TY_B21*TY_G113*TY_G29 - TY_B25*TY_G19*TY_G29 + 2*TY_B14*TY_G210*TY_G29 + 2*TY_B12*TY_G212*TY_G29
1137       
1138        TY_w[16] = -(TY_B24*TY_G110*TY_G210) - TY_B23*TY_G111*TY_G210 - TY_B22*TY_G112*TY_G210 - TY_B21*TY_G113*TY_G210 - TY_B25*TY_G19*TY_G210 - TY_B23*TY_G110*TY_G211 - TY_B22*TY_G111*TY_G211 - TY_B21*TY_G112*TY_G211
1139        TY_w[16] += -TY_B25*TY_G18*TY_G211 - TY_B24*TY_G19*TY_G211 - TY_B22*TY_G110*TY_G212 - TY_B21*TY_G111*TY_G212 - TY_B25*TY_G17*TY_G212 - TY_B24*TY_G18*TY_G212 - TY_B23*TY_G19*TY_G212 + 2*TY_B12*TY_G210*TY_G212
1140        TY_w[16] += -TY_B21*TY_G110*TY_G213 - TY_B25*TY_G16*TY_G213 - TY_B24*TY_G17*TY_G213 - TY_B23*TY_G18*TY_G213 - TY_B22*TY_G19*TY_G213 - TY_B25*TY_G113*TY_G26 - TY_B25*TY_G112*TY_G27 - TY_B24*TY_G113*TY_G27
1141        TY_w[16] += 2*TY_B14*TY_G213*TY_G27 - TY_B25*TY_G111*TY_G28 - TY_B24*TY_G112*TY_G28 - TY_B23*TY_G113*TY_G28 + 2*TY_B14*TY_G212*TY_G28
1142        TY_w[16] += -TY_G214*(TY_B25*TY_G15 + TY_B24*TY_G16 + TY_B23*TY_G17 + TY_B22*TY_G18 + TY_B21*TY_G19 - 2*TY_B14*TY_G26 - 2*TY_B12*TY_G28) - TY_B25*TY_G110*TY_G29 - TY_B24*TY_G111*TY_G29 - TY_B23*TY_G112*TY_G29
1143        TY_w[16] += -TY_B22*TY_G113*TY_G29 + 2*TY_B14*TY_G211*TY_G29 + 2*TY_B12*TY_G213*TY_G29 + TY_B34*(2*(TY_G113*TY_G17 + TY_G112*TY_G18 + TY_G111*TY_G19) + pow(TY_G110,2))
1144        TY_w[16] += TY_B32*(2*TY_G110*TY_G112 + 2*TY_G113*TY_G19 + pow(TY_G111,2)) + TY_B14*pow(TY_G210,2) + TY_B12*pow(TY_G211,2)
1145       
1146        TY_w[17] = 2*TY_B32*(TY_G111*TY_G112 + TY_G110*TY_G113) + 2*TY_B34*(TY_G110*TY_G111 + TY_G113*TY_G18 + TY_G112*TY_G19) - TY_B25*TY_G110*TY_G210 - TY_B24*TY_G111*TY_G210 - TY_B23*TY_G112*TY_G210 - TY_B22*TY_G113*TY_G210 
1147        TY_w[17] += -TY_B24*TY_G110*TY_G211 - TY_B23*TY_G111*TY_G211 - TY_B22*TY_G112*TY_G211 - TY_B21*TY_G113*TY_G211 - TY_B25*TY_G19*TY_G211 + 2*TY_B14*TY_G210*TY_G211 - TY_B23*TY_G110*TY_G212 - TY_B22*TY_G111*TY_G212 
1148        TY_w[17] += -TY_B21*TY_G112*TY_G212 - TY_B25*TY_G18*TY_G212 - TY_B24*TY_G19*TY_G212 + 2*TY_B12*TY_G211*TY_G212 - TY_B22*TY_G110*TY_G213 - TY_B21*TY_G111*TY_G213 - TY_B25*TY_G17*TY_G213 - TY_B24*TY_G18*TY_G213 
1149        TY_w[17] += -TY_B23*TY_G19*TY_G213 + 2*TY_B12*TY_G210*TY_G213 - TY_B25*TY_G113*TY_G27 - TY_B25*TY_G112*TY_G28 - TY_B24*TY_G113*TY_G28 + 2*TY_B14*TY_G213*TY_G28 - TY_B25*TY_G111*TY_G29 - TY_B24*TY_G112*TY_G29
1150        TY_w[17] += -TY_B23*TY_G113*TY_G29 + 2*TY_B14*TY_G212*TY_G29 - TY_G214*(TY_B21*TY_G110 + TY_B25*TY_G16 + TY_B24*TY_G17 + TY_B23*TY_G18 + TY_B22*TY_G19 - 2*TY_B14*TY_G27 - 2*TY_B12*TY_G29)
1151       
1152        TY_w[18] = -(TY_B25*TY_G111*TY_G210) - TY_B24*TY_G112*TY_G210 - TY_B23*TY_G113*TY_G210 - TY_B25*TY_G110*TY_G211 - TY_B24*TY_G111*TY_G211 - TY_B23*TY_G112*TY_G211 - TY_B22*TY_G113*TY_G211 - TY_B24*TY_G110*TY_G212
1153        TY_w[18] += -TY_B23*TY_G111*TY_G212 - TY_B22*TY_G112*TY_G212 - TY_B21*TY_G113*TY_G212 - TY_B25*TY_G19*TY_G212 + 2*TY_B14*TY_G210*TY_G212 - TY_B23*TY_G110*TY_G213 - TY_B22*TY_G111*TY_G213 - TY_B21*TY_G112*TY_G213
1154        TY_w[18] += -TY_B25*TY_G18*TY_G213 - TY_B24*TY_G19*TY_G213 + 2*TY_B12*TY_G211*TY_G213 - TY_B25*TY_G113*TY_G28
1155        TY_w[18] += -TY_G214*(TY_B22*TY_G110 + TY_B21*TY_G111 + TY_B25*TY_G17 + TY_B24*TY_G18 + TY_B23*TY_G19 - 2*TY_B12*TY_G210 - 2*TY_B14*TY_G28) - TY_B25*TY_G112*TY_G29 - TY_B24*TY_G113*TY_G29 + 2*TY_B14*TY_G213*TY_G29
1156        TY_w[18] += TY_B34*(2*TY_G110*TY_G112 + 2*TY_G113*TY_G19 + pow(TY_G111,2)) + TY_B32*(2*TY_G111*TY_G113 + pow(TY_G112,2)) + TY_B14*pow(TY_G211,2) + TY_B12*pow(TY_G212,2)
1157       
1158        TY_w[19] = 2*TY_B32*TY_G112*TY_G113 + 2*TY_B34*(TY_G111*TY_G112 + TY_G110*TY_G113) - TY_B25*TY_G112*TY_G210 - TY_B24*TY_G113*TY_G210 - TY_B25*TY_G111*TY_G211 - TY_B24*TY_G112*TY_G211 - TY_B23*TY_G113*TY_G211
1159        TY_w[19] += -TY_B25*TY_G110*TY_G212 - TY_B24*TY_G111*TY_G212 - TY_B23*TY_G112*TY_G212 - TY_B22*TY_G113*TY_G212 + 2*TY_B14*TY_G211*TY_G212 - TY_B24*TY_G110*TY_G213 - TY_B23*TY_G111*TY_G213 - TY_B22*TY_G112*TY_G213
1160        TY_w[19] += -TY_B21*TY_G113*TY_G213 - TY_B25*TY_G19*TY_G213 + 2*TY_B14*TY_G210*TY_G213 + 2*TY_B12*TY_G212*TY_G213 - TY_B25*TY_G113*TY_G29
1161        TY_w[19] += -TY_G214*(TY_B23*TY_G110 + TY_B22*TY_G111 + TY_B21*TY_G112 + TY_B25*TY_G18 + TY_B24*TY_G19 - 2*TY_B12*TY_G211 - 2*TY_B14*TY_G29)
1162       
1163        TY_w[20] = -(TY_B25*TY_G113*TY_G210) - TY_B25*TY_G112*TY_G211 - TY_B24*TY_G113*TY_G211 - TY_B25*TY_G111*TY_G212 - TY_B24*TY_G112*TY_G212 - TY_B23*TY_G113*TY_G212 - TY_B25*TY_G110*TY_G213 - TY_B24*TY_G111*TY_G213
1164        TY_w[20] += -TY_B23*TY_G112*TY_G213 - TY_B22*TY_G113*TY_G213 + 2*TY_B14*TY_G211*TY_G213 - (TY_B24*TY_G110 + TY_B23*TY_G111 + TY_B22*TY_G112 + TY_B21*TY_G113 + TY_B25*TY_G19 - 2*TY_B14*TY_G210 - 2*TY_B12*TY_G212)*TY_G214
1165        TY_w[20] += TY_B34*(2*TY_G111*TY_G113 + pow(TY_G112,2)) + TY_B32*pow(TY_G113,2) + TY_B14*pow(TY_G212,2) + TY_B12*pow(TY_G213,2)
1166       
1167        TY_w[21] = TY_B25*(TY_A23*TY_B14*(-3*TY_A52*TY_B24*TY_B25 + (2*TY_A43*TY_B24 + TY_A42*TY_B25)*TY_B34) + TY_B25*(TY_A22*TY_B14*(-(TY_A52*TY_B25) + TY_A43*TY_B34) + TY_A12*(4*TY_A52*TY_B24*TY_B25 - (3*TY_A43*TY_B24 + TY_A42*TY_B25)*TY_B34)))*pow(TY_B34,3)
1168       
1169        TY_w[22] = (-(TY_A23*TY_B14) + TY_A12*TY_B25)*(TY_A52*TY_B25 - TY_A43*TY_B34)*pow(TY_B25,2)*pow(TY_B34,3)
1170       
1171        if( prnt )
1172                printf "\rCoefficients of polynomial\r"
1173                variable i
1174                for ( i = 0; i < 23; i+=1 )
1175                        printf "w[%d] = %g\r", i, TY_w[i]
1176                endfor
1177                printf "\r"
1178        endif
1179       
1180end
1181
1182Function TY_capQ( d2 )
1183        Variable d2
1184       
1185        NVAR TY_B32 = root:yuk:TY_B32
1186        NVAR TY_B34 = root:yuk:TY_B34
1187
1188        return d2 * TY_B32 + pow( d2, 3 ) *  TY_B34
1189end
1190
1191Function TY_V( d2 )
1192        variable d2
1193       
1194        NVAR TY_G13 = root:yuk:TY_G13
1195        NVAR TY_G14 = root:yuk:TY_G14
1196        NVAR TY_G15 = root:yuk:TY_G15
1197        NVAR TY_G16 = root:yuk:TY_G16
1198        NVAR TY_G17 = root:yuk:TY_G17
1199        NVAR TY_G18 = root:yuk:TY_G18
1200        NVAR TY_G19 = root:yuk:TY_G19
1201        NVAR TY_G110 = root:yuk:TY_G110
1202        NVAR TY_G111 = root:yuk:TY_G111
1203        NVAR TY_G112 = root:yuk:TY_G112
1204        NVAR TY_G113 = root:yuk:TY_G113
1205
1206        return  -( pow( d2, 2 ) * TY_G13 + pow( d2, 3 ) * TY_G14 + pow( d2, 4 ) * TY_G15 + pow( d2, 5 ) * TY_G16 + pow( d2, 6 ) * TY_G17 + pow( d2, 7 ) *  TY_G18 + pow( d2, 8 ) * TY_G19 + pow( d2, 9 ) * TY_G110 +  pow( d2, 10 ) *  TY_G111 + pow( d2, 11 ) *  TY_G112 + pow( d2, 12 ) *  TY_G113 )
1207end
1208
1209Function TY_capW( d2 )
1210        Variable d2
1211       
1212        variable tmp
1213
1214        NVAR TY_G22 = root:yuk:TY_G22
1215        NVAR TY_G23 = root:yuk:TY_G23
1216        NVAR TY_G24 = root:yuk:TY_G24
1217        NVAR TY_G25 = root:yuk:TY_G25
1218        NVAR TY_G26 = root:yuk:TY_G26
1219        NVAR TY_G27 = root:yuk:TY_G27
1220        NVAR TY_G28 = root:yuk:TY_G28
1221        NVAR TY_G29 = root:yuk:TY_G29
1222        NVAR TY_G210 = root:yuk:TY_G210
1223        NVAR TY_G211 = root:yuk:TY_G211
1224        NVAR TY_G212 = root:yuk:TY_G212
1225        NVAR TY_G213 = root:yuk:TY_G213
1226        NVAR TY_G214 = root:yuk:TY_G214
1227
1228       
1229        tmp = d2  * TY_G22 + pow( d2, 2 ) * TY_G23 + pow( d2, 3 ) * TY_G24 + pow( d2, 4 ) * TY_G25 + pow( d2, 5 ) * TY_G26
1230        tmp += pow( d2, 6 ) * TY_G27 + pow( d2, 7 ) *  TY_G28 + pow( d2, 8 ) *  TY_G29 + pow( d2, 9 ) * TY_G210
1231        tmp += pow( d2, 10 ) * TY_G211 + pow( d2, 11 ) * TY_G212 + pow( d2, 12 ) * TY_G213 + pow( d2, 13 ) * TY_G214
1232       
1233        return tmp
1234end
1235
1236Function TY_X( d2 )
1237        Variable d2
1238
1239        return TY_V( d2 ) / TY_capW( d2 )
1240end
1241
1242// solve the linear system depending on d1, d2 using Cramer's rule
1243//
1244// a,b,c1,c2 are  passed by reference and returned
1245//
1246Function TY_SolveLinearEquations(  d1,  d2, a, b, c1, c2)
1247        Variable   d1,  d2, &a, &b, &c1, &c2
1248
1249
1250        NVAR TY_q22 = root:yuk:TY_q22
1251        NVAR TY_qa12 = root:yuk:TY_qa12
1252        NVAR TY_qa21 = root:yuk:TY_qa21
1253        NVAR TY_qa22 = root:yuk:TY_qa22
1254        NVAR TY_qa23 = root:yuk:TY_qa23
1255        NVAR TY_qa32 = root:yuk:TY_qa32
1256
1257        NVAR TY_qb12 = root:yuk:TY_qb12
1258        NVAR TY_qb21 = root:yuk:TY_qb21
1259        NVAR TY_qb22 = root:yuk:TY_qb22
1260        NVAR TY_qb23 = root:yuk:TY_qb23
1261        NVAR TY_qb32 = root:yuk:TY_qb32
1262
1263        NVAR TY_qc112 = root:yuk:TY_qc112
1264        NVAR TY_qc121 = root:yuk:TY_qc121
1265        NVAR TY_qc122 = root:yuk:TY_qc122
1266        NVAR TY_qc123 = root:yuk:TY_qc123
1267        NVAR TY_qc132 = root:yuk:TY_qc132
1268
1269        NVAR TY_qc212 = root:yuk:TY_qc212
1270        NVAR TY_qc221 = root:yuk:TY_qc221
1271        NVAR TY_qc222 = root:yuk:TY_qc222
1272        NVAR TY_qc223 = root:yuk:TY_qc223
1273        NVAR TY_qc232 = root:yuk:TY_qc232
1274
1275        Variable det    = TY_q22 * d1 * d2
1276        Variable det_a  = TY_qa12  * d2 + TY_qa21  * d1 + TY_qa22  * d1 * d2 + TY_qa23  * d1 * pow( d2, 2 ) + TY_qa32  * pow( d1, 2 ) * d2
1277        Variable det_b  = TY_qb12  * d2 + TY_qb21  * d1 + TY_qb22  * d1 * d2 + TY_qb23  * d1 * pow( d2, 2 ) + TY_qb32  * pow( d1, 2 ) * d2
1278        Variable det_c1 = TY_qc112 * d2 + TY_qc121 * d1 + TY_qc122 * d1 * d2 + TY_qc123 * d1 * pow( d2, 2 ) + TY_qc132 * pow( d1, 2 ) * d2
1279        Variable det_c2 = TY_qc212 * d2 + TY_qc221 * d1 + TY_qc222 * d1 * d2 + TY_qc223 * d1 * pow( d2, 2 ) + TY_qc232 * pow( d1, 2 ) * d2
1280       
1281        a  = det_a  / det
1282        b  = det_b  / det
1283        c1 = det_c1 / det
1284        c2 = det_c2 / det
1285end
1286
1287//Solve the system of linear and nonlinear equations for given Zi, Ki, phi which gives at
1288// most 22 solutions for the parameters a,b,ci,di. From the set of solutions choose the
1289// physical one and return it.
1290//
1291//
1292// a,b,c1,c2,d1,d2 are  passed by reference and returned
1293//
1294Function TY_SolveEquations(  Z1,  Z2,  K1,  K2,  phi, a,  b,  c1,  c2,  d1,  d2, prnt )
1295        Variable   Z1,  Z2,  K1,  K2,  phi, &a, &b, &c1, &c2, &d1, &d2, prnt
1296       
1297       
1298        // reduce system to a polynomial from which all solution are extracted
1299        // by doing that a lot of global background variables are set
1300        TY_ReduceNonlinearSystem( Z1, Z2, K1, K2, phi, prnt )
1301       
1302        // the two coupled non-linear eqautions were reduced to a
1303        // 22nd order polynomial, the roots are give all possible solutions
1304        // for d2, than d1 can be computed by the function X
1305       
1306        Make/O/D/N=23 real_coefficient,imag_coefficient
1307        Make/O/D/N=22 real_root,imag_root
1308       
1309        //integer degree of polynomial
1310        variable degree = 22
1311        Variable i
1312       
1313        WAVE TY_w = TY_w
1314       
1315       
1316        ////
1317        // now I need to replace this solution with FindRoots/P to get the polynomial roots
1318        ////
1319       
1320        // vector of real and imaginary coefficients in order of INCREASING powers
1321        for ( i = 0; i <= degree; i+=1 )
1322                // the global variablw TY_w was set by TY_ReduceNonlinearSystem
1323                real_coefficient[i] = TY_w[i]
1324//              imag_coefficient[i] = 0.;
1325        endfor
1326       
1327//      zrhqr(real_coefficient, degree, NR_r, NR_i);
1328       
1329        FindRoots/P=real_coefficient
1330       
1331        WAVE/C W_polyRoots = W_polyRoots
1332       
1333        for(i=0; i<degree; i+=1)
1334                real_root[i] = real(W_polyRoots[i])
1335                imag_root[i] = imag(W_polyRoots[i])
1336        endfor
1337       
1338        //end - NR solution of polynomial
1339       
1340       
1341        // show the result if in debug mode
1342        Variable x, y
1343        if ( prnt )
1344                for ( i = 0; i < degree; i+=1 )
1345                        x = real_root[i]
1346                        y = imag_root[i]
1347                        if ( chop( y ) == 0 )
1348                                printf "root(%d) = %g\r", i+1, x
1349                        else
1350                                printf "root(%d) = %g + %g i\r", i+1, x, y
1351                        endif
1352                endfor
1353                printf "\r"
1354        endif
1355       
1356       
1357       
1358        // select real roots and those satisfying Q(x) != 0 and W(x) != 0
1359        // Paper: Cluster formation in two-Yukawa Fluids, J. Chem. Phys. 122, 2005
1360        // The right set of (a, b, c1, c2, d1, d2) should have the following properties:
1361        // (1) a > 0
1362        // (2) d1, d2 are real
1363        // (3) vi/Ki > 0 <=> g(Zi) > 0
1364        // (4) if there is still more than root, calculate g(r) for each root
1365        //     and g(r) of the correct root should have the minimum average value
1366        //         inside the hardcore 
1367        Variable var_a, var_b, var_c1, var_c2, var_d1, var_d2
1368        Make/O/D/N=22 sol_a, sol_b, sol_c1, sol_c2, sol_d1, sol_d2
1369       
1370        Variable j = 0
1371        for ( i = 0; i < degree; i+=1 )
1372       
1373                x = real_root[i]
1374                y = imag_root[i]
1375                               
1376                if ( chop( y ) == 0 && TY_capW( x ) != 0 && TY_capQ( x ) != 0 )
1377               
1378                        var_d1 = TY_X( x )
1379                        var_d2 = x
1380                       
1381                        // solution of linear system for given d1, d2 to obtain a,b,ci,di
1382                        // var_a, var_b, var_c1, var_c2 passed by reference
1383                        TY_SolveLinearEquations( var_d1, var_d2, var_a, var_b, var_c1, var_c2 )
1384                       
1385                        // select physical solutions, for details check paper: "Cluster formation in
1386                        // two-Yukawa fluids", J. Chem. Phys. 122 (2005)
1387                        if ( var_a > 0 && TY_g( Z1, phi, Z1, Z2, var_a, var_b, var_c1, var_c2, var_d1, var_d2 ) > 0 && TY_g( Z2, phi, Z1, Z2, var_a, var_b, var_c1, var_c2, var_d1, var_d2 ) > 0 )
1388                                sol_a[j]  = var_a
1389                                sol_b[j]  = var_b
1390                                sol_c1[j] = var_c1
1391                                sol_c2[j] = var_c2
1392                                sol_d1[j] = var_d1
1393                                sol_d2[j] = var_d2
1394                               
1395                                if ( prnt )
1396                                        Variable eq1 = chop( TY_LinearEquation_1( Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] ) )
1397                                        Variable eq2 = chop( TY_LinearEquation_2( Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] ) )
1398                                        Variable eq3 = chop( TY_LinearEquation_3( Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] ) )
1399                                        Variable eq4 = chop( TY_LinearEquation_4( Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] ) )
1400                                        Variable eq5 = chop( TY_NonlinearEquation_1( Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] ) )
1401                                        Variable eq6 = chop( TY_NonlinearEquation_2( Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] ) )
1402                                       
1403                                        printf "solution[%d] = (%g, %g, %g, %g, %g, %g), ( eq == 0 ) = (%g, %g, %g, %g, %g, %g)\r", j, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j], eq1 , eq2, eq3, eq4, eq5, eq6
1404                                endif
1405                               
1406                                j+=1
1407                        endif           //var_a >0...
1408                endif           //chop
1409        endfor
1410        // number  remaining roots
1411        Variable n_roots = j
1412       
1413        // if there is still more than one root left, than choose the one with the minimum
1414        // average value inside the hardcore
1415        if ( n_roots > 1 )
1416               
1417                /////
1418                // it seems like this section should all be replaced in bulk with internal FFT code, rather than slow integration
1419                //
1420                // -- also, be sure to handle r=0, or the sum will always be INF
1421                ////           
1422               
1423                // the number of q values should be a power of 2
1424                // in order to speed up the FFT
1425///             int n = 1 << 14;
1426                Variable n=16384                //2^14 points
1427               
1428                // the maximum q value should be large enough
1429                // to enable a reasoble approximation of g(r)
1430                variable qmax = 16 * 10 * 2 * pi
1431                Variable q, dq = qmax / ( n - 1 )
1432               
1433                // step size for g(r)
1434                variable dr
1435               
1436                // allocate memory for pair correlation function g(r)
1437                // and structure factor S(q)
1438                Make/O/D/N=(n) sq,gr            //gr will be redimensioned!!
1439               
1440                // loop over all remaining roots
1441                Variable minVal = 1e50          //a really big number
1442                Variable selected_root = 10     
1443                Variable sumVal = 0
1444               
1445                for ( j = 0; j < n_roots; j+=1)
1446
1447                        // calculate structure factor at different q values
1448                        for ( i = 0; i < n; i+=1)
1449                       
1450                                q = dq * i
1451                                sq[i] = SqTwoYukawa( q, Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] )   
1452                               
1453                                if(i<10 && prnt)
1454                                        printf "after SqTwoYukawa: s(q) = %g\r",sq[i]
1455                                endif
1456                               
1457                        endfor
1458                       
1459                        // calculate pair correlation function for given
1460                        // structure factor, g(r) is computed at values
1461                        // r(i) = i * dr
1462
1463//                      Yuk_SqToGr( phi, dq, sq, dr, gr, n )
1464
1465
1466                        Yuk_SqToGr_FFT( phi, dq, sq, dr, gr, n )
1467       
1468                        // determine sum inside the hardcore
1469                        // 0 =< r < 1 of the pair-correlation function
1470                        sumVal = 0
1471                        for (i = 0; i < floor( 1. / dr ); i+=1 )
1472                       
1473                                sumVal += abs( gr[i] )
1474                               
1475                                if(i<10 && prnt)
1476                                        printf "g(r) in core = %g\r",abs(gr[i])
1477                                endif
1478                               
1479                        endfor
1480
1481                        if ( sumVal < minVal )
1482                                minVal = sumVal
1483                                selected_root = j
1484                        endif
1485                       
1486                        if(prnt)
1487                                printf "min = %g  sum = %g\r",minVal,sumVal
1488                        endif
1489                       
1490                endfor 
1491
1492               
1493                // physical solution was found
1494                a  = sol_a [ selected_root ]            //sol_a [ selected_root ];
1495                b  = sol_b [ selected_root ]
1496                c1 = sol_c1[ selected_root ]
1497                c2 = sol_c2[ selected_root ]
1498                d1 = sol_d1[ selected_root ]
1499                d2 = sol_d2[ selected_root ]
1500               
1501                return 1
1502       
1503        else
1504                if ( n_roots == 1 )
1505       
1506                        a  = sol_a [0]
1507                        b  = sol_b [0]
1508                        c1 = sol_c1[0]
1509                        c2 = sol_c2[0]
1510                        d1 = sol_d1[0]
1511                        d2 = sol_d2[0]
1512                       
1513                        return 1
1514                else           
1515                        // no solution was found
1516                        return 0
1517                endif
1518        endif
1519       
1520end
1521
1522
1523
1524
1525//
1526Function Yuk_SqToGr_FFT( phi, dq, sq, dr, gr, n )
1527        Variable  phi, dq
1528        WAVE sq
1529        Variable &dr
1530        WAVE gr
1531        Variable n
1532
1533        Variable npts,ii,rval,jj,qval,spread=1
1534        Variable alpha
1535
1536       
1537        WaveStats/Q sq
1538        npts = V_npnts
1539       
1540        dr = 2*pi/(npts*dq)
1541       
1542        Make/O/D/N=(npts) temp
1543       
1544        temp = p*(sq[p] - 1)
1545        alpha = npts * pow( dq, 3 ) / ( 24 * pi * pi * phi )
1546       
1547        FFT/OUT=1/DEST=W_FFT temp
1548       
1549       
1550        WAVE/C W_FFT = W_FFT
1551
1552        Redimension/N=(numpnts(W_FFT)) gr
1553       
1554        gr = 1 + alpha/p*imag(W_FFT)
1555       
1556        gr[0] = 0
1557       
1558//      Killwaves/Z temp
1559
1560        return(0)
1561
1562End
1563
1564////////////////////////end converted procedures //////////////////////////////////
Note: See TracBrowser for help on using the repository browser.