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

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

Changed 2-yukawa to have better names for the parameters, a macro to plot g(r), and to return S(q)=1000 if the calculation fails, to be more obvious.

Added a Corrlation matrix function (in Wrapper.ipf), but it's not linked anywhere yet, since I don't know yet how best to display it.

Simplified the function to list .ASC files

Added comments to the polarization panel procedures.

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