source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2010/Two_Yukawa_v40.ipf @ 931

Last change on this file since 931 was 931, checked in by srkline, 9 years ago

Two changes to existing model functions:

1) For the 2-yukawa model: Added conditions to enforce Z1>Z2. If this condition is not met, then a valid-looking (but incorrect) S(q) may be returned. These conditions were identified and satisfied within Yun's original Matalb code but were not transferred to the Igor code or XOP. These conditions are enforced within the Igor code before sending to either the XOP or local Igor code.

2) For the Teubner-Strey model: The coefficients have been re-written in terms of the two length scales and a proper absolute scaling. The Macro to calculate the length scales now also reports the amphiphilicity factor.

3) The help file entries for both of these functions have been appropriately updated.

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