Context Navigation

Yukawa_SQ_v40.ipf @ 755

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

Added new model with Yun Liu's TwoYukawa? structure factor. A one-yukawa version is also added, but only in the XOP. Note that the SQ calculation is AAO and NOT threadsafe, unlike all other fitting functions. The SQ calculation has not yet been added to any form factors.

Corrected Fractal_PolyCore to return P*S rather than just P.

Only allow a 2D simulation to be saved inn 2D VAX format if it is RAW counts, not absolute rescaled data. Otherwise the DP values are mangled on conversion to VAX integers. In the future, the data may be "unscaled" if there is an outcry.

File size: 88.7 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	//
17	//
18	// *** passing in Z values of zero can cause the XOP to crash. test for them here and send good values.
19	// -- the XOP will be modified to handle this and noted here when it is done. 0.001 seems to be OK
20	// as a low value.
21	// -- for OneYukawa, 0.1 seems to be a reasonable minimum
22	//
23	// - remember that the dimensionless Q variable is Q*diameter
24	//
25	//
26	// conversion to Igor from the c-code was not terribly painful, and very useful for debugging.
27	//
28	//
29	//
30	// as of September 2010:
31	//
32	// the one-component has not been tested at all
33	//
34	// -- the two component result nearly matches the result that Yun gets. I do need to relax the criteria for
35	// rejecting solutions, however. The XOP code rejects solutions that Yun considers "good". I guess I
36	// need all of the intermediate values (polynomial coefficients, solution vectors, etc.). Other than some of the
37	// numerical values not matching up - the output S(q) looks to be correct.
38	//
39	// -- also, for some cases, the results are VERY finicky - ususally there is a threshold value say, in Z, where
40	// going beyond that value is unstable. Here, in can be a bit random as to which values works and which do not.
41	// It must be hitting some strange zeros in the functions.
42	//
43	//
44	// TO ADD:
45	//
46	// - a mechanism for plotting the potential, so that users have a good handle on what the parameters actually mean.
47	// - as with plotting the potential, transform the K and Z into some energy scale if possible for interpretation.
48	//
49	//
50	//
51	//
52	//
53	/////////////////////////////////////////////
54
55
56
57
58	Proc PlotOneYukawa(num,qmin,qmax)
59	Variable num=200,qmin=0.001,qmax=0.5
60	Prompt num "Enter number of data points for model: "
61	Prompt qmin "Enter minimum q-value (A^-1) for model: "
62	Prompt qmax "Enter maximum q-value (A^-1) for model: "
63
64	Make/O/D/n=(num) xwave_1yuk,ywave_1yuk
65	xwave_1yuk = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
66	Make/O/D coef_1yuk = {0.1,50,-1,10}
67	make/o/t parameters_1yuk = {"volume fraction","Radius (A)","scale, K","charge, Z"}
68	Edit parameters_1yuk,coef_1yuk
69	Variable/G root:g_1yuk
70	g_1yuk := OneYukawa(coef_1yuk,ywave_1yuk,xwave_1yuk)
71	// g_1yuk := OneYukawaX(coef_1yuk,xwave_1yuk,ywave_1yuk) //be sure to have x and y in the correct order
72	Display ywave_1yuk vs xwave_1yuk
73	ModifyGraph marker=29,msize=2,mode=4
74	Label bottom "q (A\\S-1\\M)"
75	Label left "Structure Factor"
76	AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
77
78	AddModelToStrings("OneYukawa","coef_1yuk","parameters_1yuk","1yuk")
79	End
80
81	//AAO version
82	Function OneYukawa(cw,yw,xw) : FitFunc
83	Wave cw,yw,xw
84
85	if(abs(cw[3]) < 0.1)
86	cw[3] = 0.1
87	endif
88
89	#if exists("OneYukawaX")
90	OneYukawaX(cw,xw,yw)
91	#else
92	yw = 0
93	#endif
94	return(0)
95	End
96
97
98	// no igor code, return 0
99	//
100	Function fOneYukawa(w,x) : FitFunc
101	Wave w
102	Variable x
103
104	return (0)
105	End
106
107	//////////////////////////////////////////////////////////////
108	Proc PlotTwoYukawa(num,qmin,qmax)
109	Variable num=200,qmin=0.001,qmax=0.5
110	Prompt num "Enter number of data points for model: "
111	Prompt qmin "Enter minimum q-value (A^-1) for model: "
112	Prompt qmax "Enter maximum q-value (A^-1) for model: "
113
114	declare2YGlobals() //only necessary if Igor code is used. Not needed if XOP code is used.
115
116	Make/O/D/n=(num) xwave_2yuk,ywave_2yuk
117	xwave_2yuk = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
118	Make/O/D coef_2yuk = {0.2,50,6,10,-1,2}
119	make/o/t parameters_2yuk = {"volume fraction","Radius (A)","scale, K1","charge, Z1","scale, K2","charge, Z2"}
120	Edit parameters_2yuk,coef_2yuk
121	Variable/G root:g_2yuk
122	g_2yuk := TwoYukawa(coef_2yuk,ywave_2yuk,xwave_2yuk)
123	Display ywave_2yuk vs xwave_2yuk
124	ModifyGraph marker=29,msize=2,mode=4
125	Label bottom "q (A\\S-1\\M)"
126	Label left "Structure Factor"
127	AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
128
129	AddModelToStrings("TwoYukawa","coef_2yuk","parameters_2yuk","2yuk")
130
131	End
132
133
134	//AAO version
135	Function TwoYukawa(cw,yw,xw) : FitFunc
136	Wave cw,yw,xw
137
138	if(abs(cw[2]) < 0.001)
139	cw[2] = 0.001
140	endif
141	if(abs(cw[3]) < 0.001)
142	cw[3] = 0.001
143	endif
144	if(abs(cw[4]) < 0.001)
145	cw[4] = 0.001
146	endif
147	if(abs(cw[5]) < 0.001)
148	cw[5] = 0.001
149	endif
150
151
152	#if exists("TwoYukawaX")
153	TwoYukawaX(cw,xw,yw)
154	#else
155	fTwoYukawa(cw,xw,yw)
156	#endif
157	return(0)
158	End
159
160	Proc TestTheIgor2YUK()
161	//if the regular 2-yukawa procedure is already plotted
162	// -- then append it to thte graph yourself
163	Duplicate/O ywave_2yuk ywave_2yuk_Igor
164	Variable/G root:g_2yuk_Igor=0
165	g_2yuk_Igor := fTwoYukawa(coef_2yuk,xwave_2yuk,ywave_2yuk_Igor)
166	End
167
168	//
169	Function fTwoYukawa(cw,xw,yw) : FitFunc
170	Wave cw,xw,yw
171
172	Variable Z1, Z2, K1, K2, phi,radius
173	phi = cw[0]
174	radius = cw[1]
175	K1 = cw[2]
176	Z1 = cw[3]
177	K2 = cw[4]
178	Z2 = cw[5]
179
180	Variable a,b,c1,c2,d1,d2
181
182	Variable ok,check,prnt
183	prnt = 0 //print out intermediates
184
185	ok = TY_SolveEquations( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2, prnt ) // a,b,c1,c2,d1,d2 are returned
186	if(ok)
187	check = TY_CheckSolution( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
188	if(prnt)
189	printf "solution = (%g, %g, %g, %g, %g, %g) check = %d\r", a, b, c1, c2, d1, d2, check
190	endif
191
192	// if(check)
193	if(ok) //if(ok) simply takes the best solution, not necessarily one that passes TY_CheckSolution
194	yw = SqTwoYukawa(xwradius2, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2)
195	// printf("%g %g\n",q,sq)
196	endif
197	endif
198
199	return (0)
200	End
201
202
203	///////////////////// converted procedures from c-code ////////////////////////////
204
205	/// there were two functions defined as TY_q: one as TY_Q and one as TY_q. I renamed the TY_Q function as TY_capQ, and left TY_q unchanged
206
207	// function TY_W change to TY_capW, since there is a wave named TY_w
208
209
210
211
212
213	Static Function chop(x)
214	Variable x
215
216	if ( abs(x) < 1e-6 )
217	return 0
218	else
219	return x
220	endif
221
222	end
223
224	Static Function pow(a,b)
225	Variable a,b
226
227	return (a^b)
228	end
229
230	///*
231	// ==================================================================================================
232	//
233	// The two-yukawa structure factor is uniquley determined by 6 parameters a, b, c1, c2, d1, d2,
234	// which are the solution of a system of 6 equations ( 4 linear, 2 nonlinear ). The solution can
235	// constructed by the roots of a polynomial of 22nd degree. For more details see attached
236	// Mathematica notebook, where a derivation is given
237	//
238	// ==================================================================================================
239	// */
240
241	// these all may need to be declared as global variables !!
242	//
243	// - they are defined in a global scope in the c-code!
244	//
245	// - change the data folder
246	Function declare2YGlobals()
247
248	NewDataFolder/O/S root:yuk
249
250	Variable/G TY_q22
251	Variable/G TY_qa12, TY_qa21, TY_qa22, TY_qa23, TY_qa32
252	Variable/G TY_qb12, TY_qb21, TY_qb22, TY_qb23, TY_qb32
253	Variable/G TY_qc112, TY_qc121, TY_qc122, TY_qc123, TY_qc132
254	Variable/G TY_qc212, TY_qc221, TY_qc222, TY_qc223, TY_qc232
255	Variable/G TY_A12, TY_A21, TY_A22, TY_A23, TY_A32, TY_A41, TY_A42, TY_A43, TY_A52
256	Variable/G TY_B12, TY_B14, TY_B21, TY_B22, TY_B23, TY_B24, TY_B25, TY_B32, TY_B34
257	Variable/G TY_F14, TY_F16, TY_F18, TY_F23, TY_F24, TY_F25, TY_F26, TY_F27, TY_F28, TY_F29, TY_F32, TY_F33, TY_F34, TY_F35, TY_F36, TY_F37, TY_F38, TY_F39, TY_F310
258	Variable/G TY_G13, TY_G14, TY_G15, TY_G16, TY_G17, TY_G18, TY_G19, TY_G110, TY_G111, TY_G112, TY_G113, TY_G22, TY_G23, TY_G24, TY_G25, TY_G26, TY_G27, TY_G28, TY_G29, TY_G210, TY_G211, TY_G212, TY_G213, TY_G214
259
260	SetDataFolder root:
261	//this is an array, already global TY_w[23];
262
263	End
264
265
266	Function TY_sigma( s, Z1, Z2, a, b, c1, c2, d1, d2 )
267	Variable s, Z1, Z2, a, b, c1, c2, d1, d2
268
269	return -(a / 2. + b + c1 * exp( -Z1 ) + c2 * exp( -Z2 )) / s + a * pow( s, -3 ) + b * pow( s, -2 ) + ( c1 + d1 ) * pow( s + Z1, -1 ) + ( c2 + d2 ) * pow( s + Z2, -1 )
270	end
271
272	Function TY_tau( s, Z1, Z2, a, b, c1, c2 )
273	Variable s, Z1, Z2, a, b, c1, c2
274
275	return b * pow( s, -2 ) + a * ( pow( s, -3 ) + pow( s, -2 ) ) - pow( s, -1 ) * ( c1 * Z1 * exp( -Z1 ) * pow( s + Z1, -1 ) + c2 * Z2 * exp( -Z2 ) * pow( s + Z2, -1 ) )
276	end
277
278	Function TY_q( s, Z1, Z2, a, b, c1, c2, d1, d2 )
279	Variable s, Z1, Z2, a, b, c1, c2, d1, d2
280	return TY_sigma(s, Z1, Z2, a, b, c1, c2, d1, d2) - exp( -s ) * TY_tau(s, Z1, Z2, a,b, c1, c2)
281	end
282
283	Function TY_g( s, phi, Z1, Z2, a, b, c1, c2, d1, d2 )
284	Variable s, phi, Z1, Z2, a, b, c1, c2, d1, d2
285	return s * TY_tau( s, Z1, Z2, a, b, c1, c2 ) * exp( -s ) / ( 1 - 12 * phi * TY_q( s, Z1, Z2, a, b, c1, c2, d1, d2 ) )
286	end
287
288	///*
289	// ==================================================================================================
290	//
291	// Structure factor for the potential
292	//
293	// V(r) = -kB * T * ( K1 * exp[ -Z1 * (r - 1)] / r + K2 * exp[ -Z2 * (r - 1)] / r ) for r > 1
294	// V(r) = inf for r <= 1
295	//
296	// The structure factor is parametrized by (a, b, c1, c2, d1, d2)
297	// which depend on (K1, K2, Z1, Z2, phi).
298	//
299	// ==================================================================================================
300	// */
301
302	Function TY_hq( q, Z, K, v )
303	Variable q, Z, K, v
304
305	if ( q == 0)
306	return (exp(-2.Z)(v + (v(-1. + Z) - 2.KZ)exp(Z))(-(v(1. + Z)) + (v + 2.KZ(1. + Z))exp(Z))pow(K,-1)pow(Z,-4))/4.
307	else
308
309	variable t1, t2, t3, t4
310
311	t1 = ( 1. - v / ( 2. * K * Z * exp( Z ) ) ) * ( ( 1. - cos( q ) ) / ( qq ) - 1. / ( ZZ + q*q ) )
312	t2 = ( vv ( q * cos( q ) - Z * sin( q ) ) ) / ( 4. * K * ZZ q * ( ZZ + qq ) )
313	t3 = ( q * cos( q ) + Z * sin( q ) ) / ( q * ( ZZ + qq ) )
314	t4 = v / ( Z * exp( Z ) ) - vv / ( 4. K * ZZ exp( 2. * Z ) ) - K
315
316	return v / Z * t1 - t2 + t3 * t4
317	endif
318	end
319
320
321	Function TY_pc( q, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
322	Variable q, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
323
324	variable v1 = 24. * phi * K1 * exp( Z1 ) * TY_g( Z1, phi, Z1, Z2, a, b, c1, c2, d1, d2 )
325	variable v2 = 24. * phi * K2 * exp( Z2 ) * TY_g( Z2, phi, Z1, Z2, a, b, c1, c2, d1, d2 )
326
327	variable a0 = a * a
328	variable b0 = -12. * phi ( pow( a + b,2 ) / 2. + a ( c1 * exp( -Z1 ) + c2 * exp( -Z2 ) ) )
329
330	variable t1, t2, t3
331
332	if ( q == 0 )
333	t1 = a0 / 3.
334	t2 = b0 / 4.
335	t3 = a0 * phi / 12.
336	else
337	t1 = a0 * ( sin( q ) - q * cos( q ) ) / pow( q, 3 )
338	t2 = b0 * ( 2. * q * sin( q ) - ( q * q - 2. ) * cos( q ) - 2. ) / pow( q, 4 )
339	t3 = a0 * phi * ( ( qq - 6. ) 4. * q * sin( q ) - ( pow( q, 4 ) - 12. * qq + 24.) cos( q ) + 24. ) / ( 2. * pow( q, 6 ) )
340	endif
341
342	variable t4 = TY_hq( q, Z1, K1, v1 ) + TY_hq( q, Z2, K2, v2 )
343
344	return -24. * phi * ( t1 + t2 + t3 + t4 )
345	end
346
347	Function SqTwoYukawa( q, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
348	variable q, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
349
350	if ( Z1 == Z2 )
351	// one-yukawa potential
352	return 0
353	else
354	// two-yukawa potential
355	return 1. / ( 1. - TY_pc( q, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
356	endif
357	end
358
359	///*
360	//==================================================================================================
361	//
362	// Non-linear eqaution system that determines the parameter for structure factor
363	//
364	//==================================================================================================
365	//*/
366
367	Function TY_LinearEquation_1( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
368	Variable Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
369
370	return b - 12. * phi * ( -a / 8. - b / 6. + d1 * pow( Z1, -2 ) + c1 * ( pow( Z1, -2 ) - exp( -Z1 ) * ( 0.5 + ( 1. + Z1 ) * pow( Z1, -2 ) ) ) + d2 * pow( Z2, -2 ) + c2 * ( pow( Z2, -2 ) - exp( -Z2 )* ( 0.5 + ( 1. + Z2 ) * pow( Z2, -2 ) ) ) )
371	end
372
373	Function TY_LinearEquation_2( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
374	Variable Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
375
376	return 1. - a - 12. * phi * ( -a / 3. - b / 2. + d1 * pow( Z1, -1 ) + c1 * ( pow( Z1, -1 ) - ( 1. + Z1 ) * exp( -Z1 ) * pow( Z1, -1 ) ) + d2 * pow( Z2, -1 ) + c2 * ( pow( Z2, -1 ) - ( 1. + Z2 ) * exp( -Z2 ) * pow( Z2, -1 ) ) )
377	end
378
379	Function TY_LinearEquation_3( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
380	Variable Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
381
382	return K1 * exp( Z1 ) - d1 * Z1 * ( 1. - 12. * phi * TY_q( Z1, Z1, Z2, a, b, c1, c2, d1, d2 ) )
383	end
384
385	Function TY_LinearEquation_4( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
386	Variable Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
387
388	return K2 * exp( Z2 ) - d2 * Z2 * ( 1. - 12. * phi * TY_q( Z2, Z1, Z2, a, b, c1, c2, d1, d2 ) )
389	end
390
391	Function TY_NonlinearEquation_1( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
392	Variable Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
393
394	return c1 + d1 - 12. * phi * ( ( c1 + d1 ) * TY_sigma( Z1, Z1, Z2, a, b, c1, c2, d1, d2 ) - c1 * TY_tau( Z1, Z1, Z2, a, b, c1, c2 ) * exp( -Z1 ) )
395	end
396
397	Function TY_NonlinearEquation_2( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
398	Variable Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
399
400	return c2 + d2 - 12. * phi * ( ( c2 + d2 ) * TY_sigma( Z2, Z1, Z2, a, b, c1, c2, d1, d2 ) - c2 * TY_tau( Z2, Z1, Z2, a, b, c1, c2 ) * exp( -Z2 ) )
401	end
402
403	// Check the computed solutions satisfy the system of equations
404	Function TY_CheckSolution( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
405	variable Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
406
407	variable eq_1 = chop( TY_LinearEquation_1( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
408	variable eq_2 = chop( TY_LinearEquation_2( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
409	variable eq_3 = chop( TY_LinearEquation_3( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
410	variable eq_4 = chop( TY_LinearEquation_4( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
411	variable eq_5 = chop( TY_NonlinearEquation_1( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
412	variable eq_6 = chop( TY_NonlinearEquation_2( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
413
414	// printf("Check of solution = %g %g %g %g %g %g\r",eq_1,eq_2,eq_3,eq_4,eq_5,eq_6);
415	// check if all equation are zero
416	return ( eq_1 == 0 && eq_2 == 0 && eq_3 == 0 && eq_4 == 0 && eq_5 == 0 && eq_6 == 0 )
417	end
418
419	Function TY_ReduceNonlinearSystem( Z1, Z2, K1, K2, phi, prnt )
420	Variable Z1, Z2, K1, K2, phi, prnt
421
422
423	// /* solution of the 4 linear equations depending on d1 and d2, the solution is polynomial
424	// in d1, d2. We represend the solution as determiants obtained by Cramer's rule
425	// which can be expressed by their coefficient matrices
426	// */
427
428	Variable m11 = (3.*phi)/2.
429	Variable m13 = 6.phiexp(-Z1)(2. + Z1(2. + Z1) - 2.exp(Z1))pow(Z1,-2)
430	Variable m14 = 6.phiexp(-Z2)(2. + Z2(2. + Z2) - 2.exp(Z2))pow(Z2,-2)
431	Variable m23 = -12.phiexp(-Z1)(-1. - Z1 + exp(Z1))pow(Z1,-1)
432	Variable m24 = -12.phiexp(-Z2)(-1. - Z2 + exp(Z2))pow(Z2,-1)
433	Variable m31 = -6.phiexp(-Z1)pow(Z1,-2)(2.(1 + Z1) + exp(Z1)(-2. + pow(Z1,2)))
434	Variable m32 = -12.phi(-1. + Z1 + exp(-Z1))*pow(Z1,-1)
435	Variable m33 = 6.phiexp(-2.Z1)pow(-1. + exp(Z1),2)
436	Variable m34 = 12.phiexp(-Z1 - Z2)(Z2 - (Z1 + Z2)exp(Z1) + Z1exp(Z1 + Z2))pow(Z1 + Z2,-1)
437	Variable m41 = -6.phiexp(-Z2)pow(Z2,-2)(2.(1. + Z2) + exp(Z2)(-2. + pow(Z2,2)))
438	Variable m42 = -12.phi(-1. + Z2 + exp(-Z2))*pow(Z2,-1)
439	Variable m43 = 12.phiexp(-Z1 - Z2)(Z1 - (Z1 + Z2 - Z2exp(Z1))exp(Z2))pow(Z1 + Z2,-1)
440	Variable m44 = 6.phiexp(-2Z2)pow(-1. + exp(Z2),2)
441
442	// /* determinant of the linear system expressed as coefficient matrix in d1, d2 */
443
444	NVAR TY_q22 = root:yuk:TY_q22
445
446	TY_q22 = m14(-(m33m42) + m23(m32m41 - m31m42) + m32m43 + (4.m11(-3.m33m41 + 2.m33m42 + 3.m31m43 - 2.m32m43))/3.)
447	TY_q22 += m13(m34m42 + m24(-(m32m41) + m31m42) - m32m44 + (4.m11(3.m34m41 - 2.m34m42 - 3.m31m44 + 2.m32m44))/3.)
448	TY_q22 += (3.m24(m33(3.m41 + 4.m11m41 - 3.m11m42) + (-3.m31 - 4.m11m31 + 3.m11m32)m43) + 3.m23(-3.m34m41 - 4.m11m34m41 + 3.m11m34m42 + 3.m31m44 + 4.m11m31m44 - 3.m11m32m44) - (m34m43 - m33m44)pow(3. - 2.m11,2))/9.
449
450	if( prnt )
451	printf "\rDet = \r"
452	// printf "%f\t%f\r%f\t%f\r", 0., 0., 0., TY_q22
453	printf "TY_q22 = %15.12g\r",TY_q22
454	endif
455
456	// /* Matrix representation of the determinant of the of the system where row refering to
457	// the variable a is replaced by solution vector */
458
459	NVAR TY_qa12 = root:yuk:TY_qa12
460	NVAR TY_qa21 = root:yuk:TY_qa21
461	NVAR TY_qa22 = root:yuk:TY_qa22
462	NVAR TY_qa23 = root:yuk:TY_qa23
463	NVAR TY_qa32 = root:yuk:TY_qa32
464
465	Variable t1,t2,t3,t4,t5,t6,t7,t8,t9,t10
466	Variable t11,t12,t13,t14,t15,t16,t17,t18,t19,t20 //simply to keep the line length small enough
467
468	TY_qa12 = (K1(3.m14(m23m42 - 4.m11m43) - 3.m13(m24m42 - 4.m11m44) + (3. + 4.m11)(m24m43 - m23m44))exp(Z1))/3.
469
470	TY_qa21 = -(K2(3.m14(m23m32 - 4.m11m33) - 3.m13(m24m32 - 4.m11m34) + (3. + 4.m11)(m24m33 - m23m34))exp(Z2))/3.
471
472	TY_qa22 = m14(-(m23m42Z1) + 4.m11m43Z1 - m33(m42 + 4.m11Z2) + m32(m43 + m23Z2)) + (3.m13(m24m42Z1 - 4.m11m44Z1 + m34(m42 + 4.m11Z2) - m32(m44 + m24Z2)) + (3. + 4.m11)(-(m24m43Z1) + m23m44Z1 - m34(m43 + m23Z2) + m33(m44 + m24*Z2)))/3.
473
474
475	t1 = (2.(-3.m13m42 + 3.m43 + 4.m11m43)Z1pow(Z2,2) - m33(Z1 + Z2)(6.m42 + (3. + 4.m11)pow(Z2,2)) + 3.m32(Z1 + Z2)(2.m43 + m13pow(Z2,2)))
476	t2 = (2.(3.m14m42 - 3.m44 - 4.m11m44)Z1pow(Z2,2) + m34(Z1 + Z2)(6.m42 + (3. + 4.m11)pow(Z2,2)) - 3.m32(Z1 + Z2)(2.m44 + m14pow(Z2,2)))
477	t3 = (3.(m14m33m42 - m13m34m42 - m14m32m43 + m34m43 + m13m32m44 - m33m44)Z2(Z1 + Z2) + 2.m11(6.(-(m14m43) + m13m44)Z1pow(Z2,2) + m34(Z1 + Z2)(2.m43(-3. + Z2) - 3.m13pow(Z2,2)) + m33(Z1 + Z2)(6.m44 - 2.m44Z2 + 3.m14*pow(Z2,2))))
478
479	TY_qa23 = 2.phipow(Z2,-2)(m24t1 + m23t2 + 2.t3)*pow(Z1 + Z2,-1)
480
481
482
483	t1 = ((-3.m13m42 + (3. + 4.m11)m43)(Z1 + Z2)pow(Z1,2) - 2.m33(3.m42(Z1 + Z2) + (3. + 4.m11)Z2pow(Z1,2)) + 6.m32(m43(Z1 + Z2) + m13Z2pow(Z1,2)))
484	t2 = ((3.m14m42 - (3. + 4.m11)m44)(Z1 + Z2)pow(Z1,2) + m34(6.m42(Z1 + Z2) + 2.(3. + 4.m11)Z2pow(Z1,2)) - 6.m32(m44(Z1 + Z2) + m14Z2pow(Z1,2)))
485	t3 = (3.(m14m33m42 - m13m34m42 - m14m32m43 + m34m43 + m13m32m44 - m33m44)Z1(Z1 + Z2) + 2.m11(-3.(m14m43 - m13m44)(Z1 + Z2)pow(Z1,2) + 2.m34(m43(-3 + Z1)(Z1 + Z2) - 3.m13Z2pow(Z1,2)) + m33(-2.m44(-3. + Z1)(Z1 + Z2) + 6.m14Z2pow(Z1,2))))
486
487	TY_qa32 = 2.phipow(Z1,-2)(m24t1 + m23t2 + 2.t3)*pow(Z1 + Z2,-1)
488
489	if( prnt )
490	printf "\rDet_a = \r"
491	// printf "%f\t%f\t%f\r%f\t%f\t%f\r%f\t%f\t%f\r", 0., TY_qa12, 0., TY_qa21, TY_qa22, TY_qa23, 0., TY_qa32, 0.
492	printf "TY_qa12 = %15.12g\r",TY_qa12
493	printf "TY_qa21 = %15.12g\r",TY_qa21
494	printf "TY_qa22 = %15.12g\r",TY_qa22
495	printf "TY_qa23 = %15.12g\r",TY_qa23
496	printf "TY_qa32 = %15.12g\r",TY_qa32
497	endif
498
499	// /* Matrix representation of the determinant of the of the system where row refering to
500	// the variable b is replaced by solution vector */
501
502	NVAR TY_qb12 = root:yuk:TY_qb12
503	NVAR TY_qb21 = root:yuk:TY_qb21
504	NVAR TY_qb22 = root:yuk:TY_qb22
505	NVAR TY_qb23 = root:yuk:TY_qb23
506	NVAR TY_qb32 = root:yuk:TY_qb32
507
508	TY_qb12 = (K1(-3.m11m24m43 + m14(-3.m23m41 + (-3. + 8.m11)m43) + 3.m11m23m44 + m13(3.m24m41 + 3.m44 - 8.m11m44))*exp(Z1))/3.
509
510	TY_qb21 = (K2(-3.m13m24m31 + 3.m11m24m33 + m14(3.m23m31 + (3. - 8.m11)m33) - 3.m13m34 + 8.m11m13m34 - 3.m11m23m34)*exp(Z2))/3.
511
512	TY_qb22 = m13(m31m44 - m24m41Z1 - m44Z1 + (8.m11m44Z1)/3. + m24m31Z2 + m34(-m41 + Z2 - (8.m11Z2)/3.)) + m14(m23m41Z1 + m43Z1 - (8.m11m43Z1)/3. + m33(m41 - Z2 + (8.m11Z2)/3.) - m31(m43 + m23Z2)) + m11(m24m43Z1 - m23m44Z1 + m34(m43 + m23Z2) - m33(m44 + m24Z2))
513
514	t1 = (-(m14m33m41) + m13m34m41 + m14m31m43 - m11m34m43 - m13m31m44 + m11m33m44)
515	t2 = (-3.m11m24m43 + m14(-3.m23m41 + (-3. + 8.m11)m43) + 3.m11m23m44 + m13(3.m24m41 + 3.m44 - 8.m11*m44))
516	t3 = (3.m24(m33m41 - m31m43) + m23(-3.m34m41 + 3.m31m44) + (-3. + 8.m11)(m34m43 - m33*m44))
517
518	TY_qb23 = 2.phi(3.m14m23m31 - 3.m13m24m31 + 3.m14m33 - 8.m11m14m33 + 3.m11m24m33 - 3.m13m34 + 8.m11m13m34 - 3.m11m23m34 + 2.t3 pow(Z2,-2) + 6.t1pow(Z2,-1) + 2.t2Z1*pow(Z1 + Z2,-1))
519
520
521	t1 = (-(m34(m23m41 + m43)) + m24(m33m41 - m31m43) + (m23m31 + m33)*m44)
522	t2 = (-(m14m33m41) + m13m34m41 + m14m31m43 - m13m31m44)
523	t3 = (m14(2.m23m31 + 2.m33 - m23m41 - m43) + m13(-2.m34 + m24(-2.*m31 + m41) + m44))
524	t4 = (16.m34m43 - 16.m33m44 - 6.m34m43Z1 + 6.m33m44Z1 + (6.m24m33 - 3.m24m43 + 8.m14(-2.m33 + m43) + (8.m13 - 3.m23)(2.m34 - m44))pow(Z1,2))
525	t5 = (2.m34m43(8. - 3.Z1) + 2.m33m44(-8. + 3.Z1) + (8.m14m43 - 3.m24m43 - 8.m13m44 + 3.m23m44)*pow(Z1,2))
526
527	TY_qb32 = 2.phipow(Z1,-2)(6.t1 + 6.t2Z1 + 3.t3pow(Z1,2) + (m11Z2t4 + m11Z1t5)* pow(Z1 + Z2,-1) + 6.(-(m14(m23m31 + m33)) + m13(m24m31 + m34))pow(Z1,3)*pow(Z1 + Z2,-1))
528
529
530	if( prnt )
531	printf "\rDet_b = \r"
532	// printf "%f\t%f\t%f\r%f\t%f\t%f\r%f\t%f\t%f\r", 0., TY_qb12, 0., TY_qb21, TY_qb22, TY_qb23, 0., TY_qb32, 0.
533	printf "TY_qb12 = %15.12g\r",TY_qb12
534	printf "TY_qb21 = %15.12g\r",TY_qb21
535	printf "TY_qb22 = %15.12g\r",TY_qb22
536	printf "TY_qb23 = %15.12g\r",TY_qb23
537	printf "TY_qb32 = %15.12g\r",TY_qb32
538	endif
539
540	// /* Matrix representation of the determinant of the of the system where row refering to
541	// the variable c1 is replaced by solution vector */
542	NVAR TY_qc112 = root:yuk:TY_qc112
543	NVAR TY_qc121 = root:yuk:TY_qc121
544	NVAR TY_qc122 = root:yuk:TY_qc122
545	NVAR TY_qc123 = root:yuk:TY_qc123
546	NVAR TY_qc132 = root:yuk:TY_qc132
547
548	TY_qc112 = -(K1exp(Z1)(9.m24m41 - 9.m14m42 + 3.m11(-12.m14m41 + 4.m24m41 + 8.m14m42 - 3.m24m42) + m44pow(3. - 2.m11,2)))/9.
549
550	TY_qc121 = (K2exp(Z2)(9.m24m31 - 9.m14m32 + 3.m11(-12.m14m31 + 4.m24m31 + 8.m14m32 - 3.m24m32) + m34pow(3. - 2.m11,2)))/9.
551
552	TY_qc122 = m14(-4.m11m41Z1 - m42Z1 + (8.m11m42Z1)/3. + m32(-m41 + Z2 - (8.m11Z2)/3.) + m31(m42 + 4.m11Z2)) + (3.m34((3. + 4.m11)m41 - 3.m11m42) + 9.m11m32m44 + 9.m24m41Z1 + 12.m11m24m41Z1 - 9.m11m24m42Z1 + 9.m44Z1 - 12.m11m44Z1 + 9.m11m24m32Z2 - 3.(3. + 4.m11)m31(m44 + m24Z2) - m34Z2pow(3. - 2.m11,2) + 4.m44Z1pow(m11,2))/9.
553
554
555	t1 = (m34(Z1 + Z2)(2.m42 + Z2(-2.m41 + Z2)) - m32(Z1 + Z2)(2.m44 + m14Z2(-2.m41 + Z2)) - 2.(m14m42 - m44)Z2(-(Z1Z2) + m31*(Z1 + Z2)))
556	t2 = (2.(3.m41 + 4.m11m41 - 3.m11m42)Z1pow(Z2,2) + 3.m32(Z1 + Z2)(2.m41 + m11pow(Z2,2)) - m31(Z1 + Z2)(6.m42 + (3. + 4.m11)pow(Z2,2)))
557	t3 = (8.m42 + 4.m41(-3. + Z2) - 3.m42Z2 + 2.pow(Z2,2))
558	t4 = (6.m44 - 2.m44Z2 + 3.m14*pow(Z2,2))
559	t5 = (-8.m32m44Z1 + m32m44(-8. + 3.Z1)Z2 + (3.m32m44 - 4.(m14(m32 + 3.m41 - 2.m42) + m44)Z1)pow(Z2,2) + m34(Z1 + Z2)t3 + 2.m31(Z1 + Z2)t4 - 4.m14m32*pow(Z2,3))
560
561	TY_qc123 = (2.phipow(Z2,-2)(9.t1 + 4.(-2.m44Z1 + m34(Z1 + Z2))pow(m11,2)pow(Z2,2) - 3.m24t2 - 6.m11t5)*pow(Z1 + Z2,-1))/3.;
562
563
564	t1 = ((m14m42 - m44)(2.m31 - Z1)Z1(Z1 + Z2) - 2.m34(m42(Z1 + Z2) - Z1(-(Z1Z2) + m41(Z1 + Z2))) + 2.m32(m44(Z1 + Z2) - m14Z1(-(Z1Z2) + m41(Z1 + Z2))))
565	t2 = (((3. + 4.m11)m41 - 3.m11m42)(Z1 + Z2)pow(Z1,2) + 6.m32(m41(Z1 + Z2) + m11Z2pow(Z1,2)) - 2.m31(3.m42(Z1 + Z2) + (3. + 4.m11)Z2pow(Z1,2)))
566	t3 = (-8.m32m44 + m34(m42(8. - 3.Z1) + 4.m41(-3. + Z1)) - 4.m31m44(-3. + Z1) + 3.m32m44Z1 - 2.(3.m14m41 - 2.m14m42 + m44)*pow(Z1,2))
567	t4 = (4.(3.m31 - 2.m32)m44 + Z1(-4.m31m44 + 3.m32m44 - 2.(m14(-6.m31 + 4.m32 + 3.m41 - 2.m42) + m44)Z1) + m34(m42(8. - 3.Z1) + 4.m41(-3. + Z1) + 4.pow(Z1,2)))
568
569	TY_qc132 = (-2.phipow(Z1,-2)(9.t1 + 4.(-2.m34Z2 + m44(Z1 + Z2))pow(m11,2)pow(Z1,2) + 3.m24t2 + 6.m11(Z1t3 + Z2t4))*pow(Z1 + Z2,-1))/3.;
570
571
572	if( prnt )
573	printf "\rDet_c1 = \r"
574	// printf "%f\t%f\t%f\r%f\t%f\t%f\r%f\t%f\t%f\r", 0., TY_qc112, 0., TY_qc121, TY_qc122, TY_qc123, 0., TY_qc132, 0.
575	printf "TY_qc112 = %15.12g\r",TY_qc112
576	printf "TY_qc121 = %15.12g\r",TY_qc121
577	printf "TY_qc122 = %15.12g\r",TY_qc122
578	printf "TY_qc123 = %15.12g\r",TY_qc123
579	printf "TY_qc132 = %15.12g\r",TY_qc132
580	endif
581
582	// /* Matrix representation of the determinant of the of the system where row refering to
583	// the variable c1 is replaced by solution vector */
584	NVAR TY_qc212 = root:yuk:TY_qc212
585	NVAR TY_qc221 = root:yuk:TY_qc221
586	NVAR TY_qc222 = root:yuk:TY_qc222
587	NVAR TY_qc223 = root:yuk:TY_qc223
588	NVAR TY_qc232 = root:yuk:TY_qc232
589
590	TY_qc212 = (K1exp(Z1)(9m23m41 - 9m13m42 + 3m11(-12m13m41 + 4m23m41 + 8m13m42 - 3m23m42) + m43pow(3 - 2m11,2)))/9.
591
592	TY_qc221 = -(K2exp(Z2)(9m23m31 - 9m13m32 + 3m11(-12m13m31 + 4m23m31 + 8m13m32 - 3m23m32) + m33pow(3 - 2m11,2)))/9.
593
594	TY_qc222 = m13(4m11m41Z1 + m42Z1 - (8m11m42Z1)/3. + m32(m41 - Z2 + (8m11Z2)/3.) - m31(m42 + 4m11Z2)) + (9m31m43 - 9(m23m41 + m43)Z1 + 9m23m31Z2 + 3m11((-4m23m41 + 3m23m42 + 4m43)Z1 + 4m31(m43 + m23Z2) - 3m32(m43 + m23Z2)) + m33(-3(3 + 4m11)m41 + 9m11m42 + Z2pow(3 - 2m11,2)) - 4m43Z1*pow(m11,2))/9.
595
596
597	t1 = (-(m33(Z1 + Z2)(2m42 + Z2(-2m41 + Z2))) + m32(Z1 + Z2)(2m43 + m13Z2(-2m41 + Z2)) + 2(m13m42 - m43)Z2(-(Z1Z2) + m31*(Z1 + Z2)))
598	t2 = (2(3m41 + 4m11m41 - 3m11m42)Z1pow(Z2,2) + 3m32(Z1 + Z2)(2m41 + m11pow(Z2,2)) - m31(Z1 + Z2)(6m42 + (3 + 4m11)pow(Z2,2)))
599	t3 = (-8m32m43Z1 + m32m43(-8 + 3Z1)Z2 + (3m32m43 - 4(m13(m32 + 3m41 - 2m42) + m43)Z1)pow(Z2,2) + m33(Z1 + Z2)(8m42 + 4m41(-3 + Z2) - 3m42Z2 + 2pow(Z2,2)) + 2m31(Z1 + Z2)(6m43 - 2m43Z2 + 3m13pow(Z2,2)) - 4m13m32pow(Z2,3))
600
601	TY_qc223 = (2phipow(Z2,-2)(9t1 - 4(-2m43Z1 + m33(Z1 + Z2))pow(m11,2)pow(Z2,2) + 3m23t2 + 6m11t3)*pow(Z1 + Z2,-1))/3.
602
603
604	t1 = ((m13m42 - m43)(2m31 - Z1)Z1(Z1 + Z2) - 2m33(m42(Z1 + Z2) - Z1(-(Z1Z2) + m41(Z1 + Z2))) + 2m32(m43(Z1 + Z2) - m13Z1(-(Z1Z2) + m41(Z1 + Z2))))
605	t2 = (((3 + 4m11)m41 - 3m11m42)(Z1 + Z2)pow(Z1,2) + 6m32(m41(Z1 + Z2) + m11Z2pow(Z1,2)) - 2m31(3m42(Z1 + Z2) + (3 + 4m11)Z2pow(Z1,2)))
606	t3 = (-8m32m43 + m33(m42(8 - 3Z1) + 4m41(-3 + Z1)) - 4m31m43(-3 + Z1) + 3m32m43Z1 - 2(3m13m41 - 2m13m42 + m43)*pow(Z1,2))
607	t4 = (4(3m31 - 2m32)m43 + Z1(-4m31m43 + 3m32m43 - 2(m13(-6m31 + 4m32 + 3m41 - 2m42) + m43)Z1) + m33(m42(8 - 3Z1) + 4m41(-3 + Z1) + 4pow(Z1,2)))
608
609	TY_qc232 = (2phipow(Z1,-2)(9t1 + 4(-2m33Z2 + m43(Z1 + Z2))pow(m11,2)pow(Z1,2) + 3m23t2 + 6m11(Z1t3 + Z2t4))*pow(Z1 + Z2,-1))/3.
610
611
612	if( prnt )
613	printf "\rDet_c2 = \r"
614	// printf "%f\t%f\t%f\r%f\t%f\t%f\r%f\t%f\t%f\r", 0., TY_qc212, 0., TY_qc221, TY_qc222, TY_qc223, 0., TY_qc232, 0.
615	printf "TY_qc212 = %15.12g\r",TY_qc212
616	printf "TY_qc221 = %15.12g\r",TY_qc221
617	printf "TY_qc222 = %15.12g\r",TY_qc222
618	printf "TY_qc223 = %15.12g\r",TY_qc223
619	printf "TY_qc232 = %15.12g\r",TY_qc232
620	endif
621
622	// /* coefficient matrices of nonlinear equation 1 */
623	NVAR TY_A12 = root:yuk:TY_A12
624	NVAR TY_A21 = root:yuk:TY_A21
625	NVAR TY_A22 = root:yuk:TY_A22
626	NVAR TY_A23 = root:yuk:TY_A23
627	NVAR TY_A32 = root:yuk:TY_A32
628	NVAR TY_A41 = root:yuk:TY_A41
629	NVAR TY_A42 = root:yuk:TY_A42
630	NVAR TY_A43 = root:yuk:TY_A43
631	NVAR TY_A52 = root:yuk:TY_A52
632
633	t1 = (Z1(2TY_qb12(-1 + Z1)(Z1 + Z2) - Z1(2TY_qc212Z1 + TY_qc112(Z1 + Z2))) + TY_qa12(Z1 + Z2)(-2 + pow(Z1,2)))
634	t2 = (exp(2Z1)t1 - TY_qc112(Z1 + Z2)pow(Z1,2) + 2(Z1 + Z2)exp(Z1)(TY_qa12 + (TY_qa12 + TY_qb12)Z1 + TY_qc112*pow(Z1,2)))
635
636	TY_A12 = 6phiTY_qc112exp(-2Z1 - Z2)pow(TY_q22,-2)pow(Z1,-3)(2TY_qc212exp(Z1)(-Z2 + (Z1 + Z2)exp(Z1))pow(Z1,2) + exp(Z2)t2)pow(Z1 + Z2,-1);
637
638
639	t1 = (2Z1(TY_qb21TY_qc112(-1 + Z1)(Z1 + Z2) + TY_qb12TY_qc121(-1 + Z1)(Z1 + Z2) - Z1(TY_qc121TY_qc212Z1 + TY_qc112(TY_qc121 + TY_qc221)Z1 + TY_qc112TY_qc121Z2)) + TY_qa21TY_qc112(Z1 + Z2)(-2 + pow(Z1,2)) + TY_qa12TY_qc121(Z1 + Z2)*(-2 + pow(Z1,2)))
640	t2 = (TY_qb21TY_qc112 + TY_qc121(TY_qa12 + TY_qb12 + 2TY_qc112Z1))
641	t3 = (2(TY_qa12TY_qc121 + TY_qa21TY_qc112(1 + Z1) + Z1t2)(Z1 + Z2)exp(Z1) + exp(2Z1)t1 - 2TY_qc112TY_qc121(Z1 + Z2)*pow(Z1,2))
642
643	TY_A21 = 6phiexp(-2Z1 - Z2)pow(TY_q22,-2)pow(Z1,-3)(2(TY_qc121TY_qc212 + TY_qc112TY_qc221)exp(Z1)(-Z2 + (Z1 + Z2)exp(Z1))pow(Z1,2) + exp(Z2)t3)*pow(Z1 + Z2,-1);
644
645
646	t1 = (TY_qb22TY_qc112 + TY_qc122(TY_qa12 + TY_qb12 + 2TY_qc112Z1))
647	t2 = (2Z1(TY_qb22TY_qc112(-1 + Z1)(Z1 + Z2) + TY_qb12TY_qc122(-1 + Z1)(Z1 + Z2) - Z1(TY_qc122TY_qc212Z1 + TY_qc112(TY_qc122 + TY_qc222)Z1 + TY_qc112TY_qc122Z2)) + TY_qa22TY_qc112(Z1 + Z2)(-2 + pow(Z1,2)) + TY_qa12TY_qc122(Z1 + Z2)*(-2 + pow(Z1,2)))
648	t3 = (12phi(TY_qa12TY_qc122 + TY_qa22TY_qc112(1 + Z1) + Z1t1)(Z1 + Z2)exp(Z1) - 2phiTY_qc112TY_qc122(Z1 + Z2)pow(Z1,2) + exp(2Z1)(6phit2 + TY_q22TY_qc112(Z1 + Z2)pow(Z1,3)))
649
650	TY_A22 = exp(-2Z1 - Z2)pow(TY_q22,-2)pow(Z1,-3)(12phi(TY_qc122TY_qc212 + TY_qc112TY_qc222)exp(Z1)(-Z2 + (Z1 + Z2)exp(Z1))pow(Z1,2) + exp(Z2)t3)pow(Z1 + Z2,-1);
651
652
653	t1 = ((TY_q22TY_qc112 + TY_qc123(TY_qc112 + TY_qc212) + TY_qc112TY_qc223)Z1 + TY_qc112TY_qc123Z2)
654	t2 = (TY_qa12TY_qc123 + TY_qa23TY_qc112(1 + Z1) + Z1(TY_qb23TY_qc112 + TY_qc123(TY_qa12 + TY_qb12 + 2TY_qc112Z1)))
655	t3 = (2Z1(TY_qb23TY_qc112(-1 + Z1)(Z1 + Z2) + TY_qb12TY_qc123(-1 + Z1)(Z1 + Z2) - Z1t1) + TY_qa23TY_qc112(Z1 + Z2)(-2 + pow(Z1,2)) + TY_qa12TY_qc123(Z1 + Z2)*(-2 + pow(Z1,2)))
656
657	TY_A23 = 6phiexp(-2Z1 - Z2)pow(TY_q22,-2)pow(Z1,-3)(2(TY_qc123TY_qc212 + TY_qc112TY_qc223)exp(Z1)(-Z2 + (Z1 + Z2)exp(Z1))pow(Z1,2) + exp(Z2)(2t2(Z1 + Z2)exp(Z1) + exp(2Z1)t3 - 2TY_qc112TY_qc123(Z1 + Z2)pow(Z1,2)))pow(Z1 + Z2,-1);
658
659
660	t1 = (TY_qb32TY_qc112 + (TY_qa23 + TY_qb23)TY_qc121 + (TY_qa21 + TY_qb21)TY_qc123 + (TY_qa12 + TY_qb12)TY_qc132 + TY_q22TY_qc112Z1 + 2(TY_qc121TY_qc123 + TY_qc112TY_qc132)Z1 + TY_qc122(TY_qa22 + TY_qb22 + TY_qc122Z1))
661	t2 = (TY_qc132TY_qc212 + TY_qc123TY_qc221 + TY_qc122TY_qc222 + TY_qc121TY_qc223 + TY_qc112*TY_qc232)
662	t3 = ((TY_q22 + TY_qc132)TY_qc212 + TY_qc123TY_qc221 + TY_qc122TY_qc222 + TY_qc121TY_qc223 + TY_qc112*TY_qc232)
663	t4 = (2TY_qc121TY_qc123 + 2TY_qc112TY_qc132 + pow(TY_qc122,2))
664	t5 = (6phi(2Z1(TY_qb12(-1 + Z1)(Z1 + Z2) - Z1((TY_qc112 + TY_qc121 + TY_qc212)Z1 + TY_qc112Z2)) + TY_qa12(Z1 + Z2)(-2 + pow(Z1,2))) + TY_qc122(Z1 + Z2)*pow(Z1,3))
665	t6 = (-2(TY_qa22TY_qc122 + TY_qa21TY_qc123 + TY_qa12TY_qc132) - 2(TY_qb32TY_qc112 + TY_qb23TY_qc121 + TY_qb22TY_qc122 + TY_qb21TY_qc123 + TY_qb12TY_qc132)Z1 + (2TY_qb32TY_qc112 + 2TY_qb23TY_qc121 + (TY_qa22 + 2TY_qb22 - TY_qc122)TY_qc122 + (TY_qa21 + 2TY_qb21 - 2TY_qc121)TY_qc123 + (TY_qa12 + 2TY_qb12 - 2TY_qc112)TY_qc132)pow(Z1,2))
666	t7 = -2TY_qa22TY_qc122Z1 - 2TY_qa21TY_qc123Z1 - 2TY_qa12TY_qc132Z1 + TY_qa32TY_qc112(Z1 + Z2)(-2 + pow(Z1,2)) + TY_qa23TY_qc121(Z1 + Z2)(-2 + pow(Z1,2)) - 2TY_qb32TY_qc112pow(Z1,2) - 2TY_qb23TY_qc121pow(Z1,2) - 2TY_qb22TY_qc122pow(Z1,2) - 2TY_qb21TY_qc123pow(Z1,2) - 2TY_qb12TY_qc132pow(Z1,2)
667	t8 = Z2t6 + 2TY_qb32TY_qc112pow(Z1,3) + 2TY_qb23TY_qc121pow(Z1,3) + TY_qa22TY_qc122pow(Z1,3) + 2TY_qb22TY_qc122pow(Z1,3) + TY_qa21TY_qc123pow(Z1,3) + 2TY_qb21TY_qc123pow(Z1,3) - 2TY_qc121TY_qc123pow(Z1,3) + TY_qa12TY_qc132pow(Z1,3)
668	t9 = (t7 + t8 + 2TY_qb12TY_qc132pow(Z1,3) - 2TY_qc112TY_qc132pow(Z1,3) - 2TY_qc132TY_qc212pow(Z1,3) - 2TY_qc123TY_qc221pow(Z1,3) - 2TY_qc122TY_qc222pow(Z1,3) - 2TY_qc121TY_qc223pow(Z1,3) - 2TY_qc112TY_qc232pow(Z1,3) - pow(TY_qc122,2)pow(Z1,3))
669	t10 = (12phi(TY_qa23TY_qc121 + TY_qa22TY_qc122 + TY_qa21TY_qc123 + TY_qa12TY_qc132 + TY_qa32TY_qc112(1 + Z1) + Z1t1)(Z1 + Z2)exp(Z1 + Z2) - 12phit2Z2exp(Z1)pow(Z1,2) + 12phit3(Z1 + Z2)exp(2Z1)pow(Z1,2) - 6phi(Z1 + Z2)exp(Z2)t4pow(Z1,2) + exp(2Z1 + Z2)(TY_q22t5 + 6phit9))
670
671	TY_A32 = exp(-2Z1 - Z2)pow(TY_q22,-2)pow(Z1,-3)t10*pow(Z1 + Z2,-1);
672
673
674	t1 = ((-(TY_qc132TY_qc221) - TY_qc121TY_qc232)Z2 + ((TY_q22 + TY_qc132)TY_qc221 + TY_qc121TY_qc232)(Z1 + Z2)*exp(Z1))
675	t2 = (TY_qa21TY_qc132 + TY_qa32TY_qc121(1 + Z1) + Z1(TY_qb32TY_qc121 + (TY_qa21 + TY_qb21)TY_qc132 + TY_qc121(TY_q22 + 2TY_qc132)*Z1))
676	t3 = (-2(TY_qa32TY_qc121 + TY_qa21(TY_q22 + TY_qc132)) - 2(TY_qb32TY_qc121 + TY_qb21(TY_q22 + TY_qc132))Z1 + (TY_q22(TY_qa21 + 2TY_qb21 - 2TY_qc121) + TY_qc121(TY_qa32 + 2TY_qb32 - 2TY_qc132) + (TY_qa21 + 2TY_qb21)TY_qc132)pow(Z1,2))
677	t4 = (-2(TY_qa32TY_qc121 + TY_qa21(TY_q22 + TY_qc132)) - 2(TY_qb32TY_qc121 + TY_qb21(TY_q22 + TY_qc132))Z1 + (TY_qa32TY_qc121 + 2TY_qb32TY_qc121 + TY_qa21TY_qc132 + 2TY_qb21TY_qc132 - 2TY_qc121TY_qc132 + TY_q22(TY_qa21 + 2TY_qb21 - 2TY_qc121 - 2TY_qc221) - 2TY_qc132TY_qc221 - 2TY_qc121TY_qc232)pow(Z1,2))
678
679	TY_A41 = 6phiexp(-2Z1 - Z2)pow(TY_q22,-2)pow(Z1,-3)(2exp(Z1)t1pow(Z1,2) + exp(Z2)(2t2(Z1 + Z2)exp(Z1) - 2TY_qc121TY_qc132(Z1 + Z2)pow(Z1,2) + exp(2Z1)(Z2t3 + Z1t4)))pow(Z1 + Z2,-1);
680
681
682	t1 = (TY_qb32TY_qc122 + (TY_qa22 + TY_qb22)TY_qc132 + TY_qc122(TY_q22 + 2TY_qc132)*Z1)
683	t2 = (TY_qc132TY_qc222 + TY_qc122TY_qc232)
684	t3 = ((TY_q22 + TY_qc132)TY_qc222 + TY_qc122TY_qc232)
685	t4 = (2Z1(TY_qb32TY_qc122(-1 + Z1)(Z1 + Z2) + TY_qb22TY_qc132(-1 + Z1)(Z1 + Z2) - Z1(TY_qc132TY_qc222Z1 + TY_qc122(TY_qc132 + TY_qc232)Z1 + TY_qc122TY_qc132Z2)) + TY_qa32TY_qc122(Z1 + Z2)(-2 + pow(Z1,2)) + TY_qa22TY_qc132(Z1 + Z2)*(-2 + pow(Z1,2)))
686	t5 = (6phit4 + (Z1 + Z2)pow(TY_q22,2)pow(Z1,3) + TY_q22(6phi(2Z1(TY_qb22(-1 + Z1)(Z1 + Z2) - Z1((TY_qc122 + TY_qc222)Z1 + TY_qc122Z2)) + TY_qa22(Z1 + Z2)(-2 + pow(Z1,2))) + TY_qc132(Z1 + Z2)pow(Z1,3)))
687	t6 = (12phi(TY_qa22TY_qc132 + TY_qa32TY_qc122(1 + Z1) + Z1t1)(Z1 + Z2)exp(Z1 + Z2) - 12phit2Z2exp(Z1)pow(Z1,2) + 12phit3(Z1 + Z2)exp(2Z1)pow(Z1,2) - 12phiTY_qc122TY_qc132(Z1 + Z2)exp(Z2)pow(Z1,2) + exp(2Z1 + Z2)*t5)
688
689	TY_A42 = exp(-2Z1 - Z2)pow(TY_q22,-2)pow(Z1,-3)t6*pow(Z1 + Z2,-1);
690
691
692	t1 = ((TY_qc132TY_qc223 + TY_qc123TY_qc232)Z2 - ((TY_q22 + TY_qc132)TY_qc223 + TY_qc123TY_qc232)(Z1 + Z2)*exp(Z1))
693	t2 = (TY_qa23TY_qc132 + TY_qa32TY_qc123(1 + Z1) + Z1(TY_qb32TY_qc123 + (TY_qa23 + TY_qb23)TY_qc132 + TY_qc123(TY_q22 + 2TY_qc132)*Z1))
694	t3 = (2TY_qa32TY_qc123 + 2TY_qa23(TY_q22 + TY_qc132) + 2(TY_qb32TY_qc123 + TY_qb23(TY_q22 + TY_qc132))Z1 - (TY_q22(TY_qa23 + 2TY_qb23 - 2TY_qc123) + TY_qc123(TY_qa32 + 2TY_qb32 - 2TY_qc132) + (TY_qa23 + 2TY_qb23)TY_qc132)*pow(Z1,2))
695	t4 = (2TY_qa32TY_qc123 + 2TY_qa23(TY_q22 + TY_qc132) + 2(TY_qb32TY_qc123 + TY_qb23(TY_q22 + TY_qc132))Z1 + (-(TY_qa32TY_qc123) - (TY_qa23 + 2TY_qb23)TY_qc132 + TY_q22(-TY_qa23 + 2(-TY_qb23 + TY_qc123 + TY_qc132 + TY_qc223)) + 2(-(TY_qb32TY_qc123) + TY_qc132(TY_qc123 + TY_qc223) + TY_qc123TY_qc232) + 2pow(TY_q22,2))*pow(Z1,2))
696
697	TY_A43 = -6phiexp(-2Z1 - Z2)pow(TY_q22,-2)pow(Z1,-3)(2exp(Z1)t1pow(Z1,2) + exp(Z2)(-2t2(Z1 + Z2)exp(Z1) + 2TY_qc123TY_qc132(Z1 + Z2)pow(Z1,2) + exp(2Z1)(Z2t3 + Z1t4)))pow(Z1 + Z2,-1);
698
699
700	t1 = (TY_qc132Z2 - (TY_q22 + TY_qc132)(Z1 + Z2)*exp(Z1))
701	t2 = (Z1(-2TY_qb32(-1 + Z1)(Z1 + Z2) + Z1((TY_q22 + TY_qc132 + 2TY_qc232)Z1 + (TY_q22 + TY_qc132)Z2)) - TY_qa32(Z1 + Z2)(-2 + pow(Z1,2)))
702	t3 = ((TY_q22 + TY_qc132)exp(2Z1)t2 + (Z1 + Z2)pow(TY_qc132,2)pow(Z1,2) - 2TY_qc132(Z1 + Z2)exp(Z1)(TY_qa32 + (TY_qa32 + TY_qb32)Z1 + (TY_q22 + TY_qc132)*pow(Z1,2)))
703
704	TY_A52 = -6phiexp(-2Z1 - Z2)pow(TY_q22,-2)pow(Z1,-3)(2TY_qc232exp(Z1)t1pow(Z1,2) + exp(Z2)t3)pow(Z1 + Z2,-1);
705
706
707	// normalize A
708	// /*double norm_A = sqrt(pow(TY_A52,2)+pow(TY_A43,2)+pow(TY_A42, 2)+pow(TY_A41, 2)+pow(TY_A32, 2)+
709	// pow(TY_A23,2)+pow(TY_A22,2)+pow(TY_A21, 2)+pow(TY_A12, 2));
710	// TY_A12 /= norm_A;
711	// TY_A21 /= norm_A;
712	// TY_A22 /= norm_A;
713	// TY_A23 /= norm_A;
714	// TY_A32 /= norm_A;
715	// TY_A41 /= norm_A;
716	// TY_A42 /= norm_A;
717	// TY_A43 /= norm_A;
718	// TY_A52 /= norm_A;*/
719
720	if( prnt )
721	printf "\rNonlinear equation 1 = \r"
722	// printf "%f\t\t%f\t\t%f\r", 0., TY_A12, 0.
723	// printf "%f\t\t%f\t\t%f\r", TY_A21, TY_A22, TY_A23
724	// printf "%f\t\t%f\t\t%f\r", 0., TY_A32, 0.
725	// printf "%f\t\t%f\t\t%f\r", TY_A41, TY_A42, TY_A43
726	// printf "%f\t\t%f\t\t%f\r", 0., TY_A52, 0.
727	printf "TY_A12 = %15.12g\r",TY_A12
728	printf "TY_A21 = %15.12g\r",TY_A21
729	printf "TY_A22 = %15.12g\r",TY_A22
730	printf "TY_A23 = %15.12g\r",TY_A23
731	printf "TY_A32 = %15.12g\r",TY_A32
732	printf "TY_A41 = %15.12g\r",TY_A41
733	printf "TY_A42 = %15.12g\r",TY_A42
734	printf "TY_A43 = %15.12g\r",TY_A43
735	printf "TY_A52 = %15.12g\r",TY_A52
736	endif
737
738	// /* coefficient matrices of nonlinear equation 2 */
739	NVAR TY_B12 = root:yuk:TY_B12
740	NVAR TY_B14 = root:yuk:TY_B14
741	NVAR TY_B21 = root:yuk:TY_B21
742	NVAR TY_B22 = root:yuk:TY_B22
743	NVAR TY_B23 = root:yuk:TY_B23
744	NVAR TY_B24 = root:yuk:TY_B24
745	NVAR TY_B25 = root:yuk:TY_B25
746	NVAR TY_B32 = root:yuk:TY_B32
747	NVAR TY_B34 = root:yuk:TY_B34
748
749
750
751	t1 = (TY_qa12TY_qc221 + TY_qa21TY_qc212(1 + Z2) + Z2(TY_qb21TY_qc212 + TY_qc221(TY_qa12 + TY_qb12 + 2TY_qc212Z2)))
752	t2 = (-(TY_qc121TY_qc212) - TY_qc112TY_qc221)
753	t3 = (TY_qb21TY_qc212(-1 + Z2)(Z1 + Z2) + TY_qb12TY_qc221(-1 + Z2)(Z1 + Z2) - Z2(TY_qc212TY_qc221Z1 + TY_qc112TY_qc221Z2 + TY_qc212(TY_qc121 + TY_qc221)*Z2))
754	t4 = (exp(Z1)(2Z2t3 + TY_qa21TY_qc212(Z1 + Z2)(-2 + pow(Z2,2)) + TY_qa12TY_qc221(Z1 + Z2)(-2 + pow(Z2,2))) + 2(TY_qc121TY_qc212 + TY_qc112TY_qc221)(Z1 + Z2)pow(Z2,2))
755
756	TY_B12 = 6phiexp(-Z1 - 2Z2)pow(TY_q22,-2)pow(Z2,-3)(-2TY_qc212TY_qc221(Z1 + Z2)exp(Z1)pow(Z2,2) + 2exp(Z2)((Z1 + Z2)t1exp(Z1) + t2Z1pow(Z2,2)) + exp(2Z2)t4)pow(Z1 + Z2,-1);
757
758
759
760
761	t1 = ((Z1 + Z2)(TY_qa12TY_qc223 + TY_qa23TY_qc212(1 + Z2) + Z2(TY_qb23TY_qc212 + (TY_qa12 + TY_qb12)TY_qc223 + TY_qc212(TY_q22 + 2TY_qc223)Z2))exp(Z1) + (-(TY_qc123TY_qc212) - TY_qc112TY_qc223)Z1*pow(Z2,2))
762	t2 = (TY_qc123TY_qc212 + TY_qc112(TY_q22 + TY_qc223))
763	t3 = (TY_qa23TY_qc212 + TY_qa12(TY_q22 + TY_qc223))
764	t4 = (TY_qa23TY_qc212 + TY_qa12(TY_q22 + TY_qc223) + (TY_qb23TY_qc212 + TY_qb12(TY_q22 + TY_qc223))*Z1)
765	t5 = (-2(TY_qb23TY_qc212 + TY_qb12(TY_q22 + TY_qc223)) + (TY_q22(TY_qa12 + 2TY_qb12 - 2TY_qc212) + TY_qc212(TY_qa23 + 2TY_qb23 - 2TY_qc223) + (TY_qa12 + 2TY_qb12)TY_qc223)Z1)
766	t6 = (TY_q22(TY_qa12 + 2TY_qb12 - 2TY_qc112 - 2TY_qc212) + TY_qc212(TY_qa23 + 2TY_qb23 - 2TY_qc123 - 2TY_qc223) + (TY_qa12 + 2TY_qb12 - 2TY_qc112)*TY_qc223)
767
768	TY_B14 = 6phiexp(-Z1 - 2Z2)pow(TY_q22,-2)pow(Z2,-3)(-2TY_qc212TY_qc223(Z1 + Z2)exp(Z1)pow(Z2,2) + 2exp(Z2)t1 + exp(2Z2)(2t2(Z1 + Z2)pow(Z2,2) + exp(Z1)(-2t3Z1 - 2t4Z2 + t5pow(Z2,2) + t6pow(Z2,3))))pow(Z1 + Z2,-1);
769
770
771
772	t1 = (TY_qc221(Z1 + Z2)exp(Z1)*pow(Z2,2))
773	t2 = (exp(Z1)(Z2(2TY_qb21(-1 + Z2)(Z1 + Z2) - Z2(2TY_qc121Z2 + TY_qc221(Z1 + Z2))) + TY_qa21(Z1 + Z2)(-2 + pow(Z2,2))) + 2TY_qc121(Z1 + Z2)pow(Z2,2))
774	t3 = (-(TY_qc121Z1pow(Z2,2)) + (Z1 + Z2)exp(Z1)(TY_qa21 + (TY_qa21 + TY_qb21)Z2 + TY_qc221pow(Z2,2)))
775
776	TY_B21 = 6phiTY_qc221exp(-Z1 - 2Z2)pow(TY_q22,-2)pow(Z2,-3)(-t1 + exp(2Z2)t2 + 2exp(Z2)t3)pow(Z1 + Z2,-1);
777
778
779
780	t1 = (TY_qb22TY_qc221 + TY_qc222(TY_qa21 + TY_qb21 + 2TY_qc221Z2))
781	t2 = ((Z1 + Z2)(TY_qa21TY_qc222 + TY_qa22TY_qc221(1 + Z2) + Z2t1)exp(Z1) + (-(TY_qc122TY_qc221) - TY_qc121TY_qc222)Z1pow(Z2,2))
782	t3 = (TY_qc122TY_qc221 + TY_qc121TY_qc222)
783	t4 = (TY_qb22TY_qc221(-1 + Z2)(Z1 + Z2) + TY_qb21TY_qc222(-1 + Z2)(Z1 + Z2) - Z2(TY_qc221TY_qc222Z1 + TY_qc121TY_qc222Z2 + TY_qc221(TY_qc122 + TY_qc222)*Z2))
784	t5 = (6phi(2Z2t4 + TY_qa22TY_qc221(Z1 + Z2)(-2 + pow(Z2,2)) + TY_qa21TY_qc222(Z1 + Z2)(-2 + pow(Z2,2))) + TY_q22TY_qc221(Z1 + Z2)*pow(Z2,3))
785
786	TY_B22 = exp(-Z1 - 2Z2)pow(TY_q22,-2)pow(Z2,-3)(-12phiTY_qc221TY_qc222(Z1 + Z2)exp(Z1)pow(Z2,2) + 12phiexp(Z2)t2 + exp(2Z2)(12phit3(Z1 + Z2)pow(Z2,2) + exp(Z1)t5))*pow(Z1 + Z2,-1);
787
788
789
790
791	t1 = (2TY_qc221TY_qc223 + 2TY_qc212TY_qc232 + pow(TY_qc222,2))
792	t2 = (TY_qb32TY_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_q22TY_qc221 + 2TY_qc221TY_qc223 + 2TY_qc212*TY_qc232 + pow(TY_qc222,2)))
793	t3 = (-(TY_qc132TY_qc212) - TY_qc123TY_qc221 - TY_qc122TY_qc222 - TY_qc121TY_qc223 - TY_qc112*TY_qc232)
794	t4 = (TY_qc132TY_qc212 + TY_qc123TY_qc221 + TY_qc122TY_qc222 + TY_qc121(TY_q22 + TY_qc223) + TY_qc112*TY_qc232)
795	t5 = (TY_qa32TY_qc212 + TY_qa23TY_qc221 + TY_qa22TY_qc222 + TY_qa21(TY_q22 + TY_qc223) + TY_qa12*TY_qc232)
796	t6 = (TY_qa32TY_qc212 + TY_qa23TY_qc221 + TY_qa22TY_qc222 + TY_qa21(TY_q22 + TY_qc223) + TY_qa12TY_qc232 + (TY_qb32TY_qc212 + TY_qb23TY_qc221 + TY_qb22TY_qc222 + TY_qb21(TY_q22 + TY_qc223) + TY_qb12TY_qc232)*Z1)
797	t7 = (TY_qb32TY_qc212 + TY_qb23TY_qc221 + TY_qb22TY_qc222 + TY_qb21(TY_q22 + TY_qc223) + TY_qb12*TY_qc232)
798	t8 = (TY_q22(TY_qa21 + 2TY_qb21 - 2TY_qc221) + (TY_qa22 + 2TY_qb22 - TY_qc222)TY_qc222 + TY_qc221(TY_qa23 + 2TY_qb23 - 2TY_qc223) + TY_qa21TY_qc223 + 2TY_qb21TY_qc223 + TY_qc212(TY_qa32 + 2TY_qb32 - 2TY_qc232) + TY_qa12TY_qc232 + 2TY_qb12*TY_qc232)
799	t9 = (TY_qa21 + 2TY_qb21 - 2TY_qc121 - 2TY_qc212 - 2TY_qc221)
800	t10 = (TY_qa22 + 2TY_qb22 - 2TY_qc122 - TY_qc222)
801	t11 = (TY_qa23 + 2TY_qb23 - 2TY_qc123 - 2*TY_qc223)
802	t12 = (TY_qa32 + 2TY_qb32 - 2TY_qc132 - 2*TY_qc232)
803	t13 = ((Z1 + Z2)exp(Z1)(TY_qa23TY_qc221 + TY_qa22TY_qc222 + TY_qa21TY_qc223 + TY_qa12TY_qc232 + TY_qa32TY_qc212(1 + Z2) + Z2t2) + t3Z1*pow(Z2,2))
804	t14 = (TY_q22t9 + t10TY_qc222 + TY_qc221t11 + (TY_qa21 + 2TY_qb21 - 2TY_qc121)TY_qc223 + TY_qc212t12 + (TY_qa12 + 2TY_qb12 - 2TY_qc112)TY_qc232)
805	t15 = (-6phi(Z1 + Z2)exp(Z1)t1pow(Z2,2) + 12phiexp(Z2)t13 + exp(2Z2)(12phit4(Z1 + Z2)pow(Z2,2) + exp(Z1)(TY_q22TY_qc222(Z1 + Z2)pow(Z2,3) - 6phi(2t5Z1 + 2t6Z2 - (-2t7 + t8Z1)pow(Z2,2) - t14pow(Z2,3)))))
806
807	TY_B23 = exp(-Z1 - 2Z2)pow(TY_q22,-2)pow(Z2,-3)t15*pow(Z1 + Z2,-1);
808
809
810	t1 = (TY_qa22TY_qc223 + TY_qa23TY_qc222(1 + Z2) + Z2(TY_qb23TY_qc222 + (TY_qa22 + TY_qb22)TY_qc223 + TY_qc222(TY_q22 + 2TY_qc223)*Z2))
811	t2 = (TY_qc123TY_qc222 + TY_qc122TY_qc223)
812	t3 = (TY_qc123TY_qc222 + TY_qc122(TY_q22 + TY_qc223))
813	t4 = (2Z2(TY_qb23TY_qc222(-1 + Z2)(Z1 + Z2) + TY_qb22TY_qc223(-1 + Z2)(Z1 + Z2) - Z2(TY_qc222TY_qc223Z1 + TY_qc122TY_qc223Z2 + TY_qc222(TY_qc123 + TY_qc223)Z2)) + TY_qa23TY_qc222(Z1 + Z2)(-2 + pow(Z2,2)) + TY_qa22TY_qc223(Z1 + Z2)*(-2 + pow(Z2,2)))
814	t5 = (6phit4 + (Z1 + Z2)pow(TY_q22,2)pow(Z2,3) + TY_q22(6phi(2Z2(TY_qb22(-1 + Z2)(Z1 + Z2) - Z2(TY_qc222Z1 + (TY_qc122 + TY_qc222)Z2)) + TY_qa22(Z1 + Z2)(-2 + pow(Z2,2))) + TY_qc223(Z1 + Z2)pow(Z2,3)))
815	t6 = (12phi(Z1 + Z2)t1exp(Z1 + Z2) - 12phiTY_qc222TY_qc223(Z1 + Z2)exp(Z1)pow(Z2,2) - 12phit2Z1exp(Z2)pow(Z2,2) + 12phit3(Z1 + Z2)exp(2Z2)pow(Z2,2) + exp(Z1 + 2Z2)*t5)
816
817	TY_B24 = exp(-Z1 - 2Z2)pow(TY_q22,-2)pow(Z2,-3)t6*pow(Z1 + Z2,-1);
818
819
820	t1 = (exp(Z1)(Z2(-2TY_qb23(-1 + Z2)(Z1 + Z2) + Z2((TY_q22 + TY_qc223)Z1 + (TY_q22 + 2TY_qc123 + TY_qc223)Z2)) - TY_qa23(Z1 + Z2)(-2 + pow(Z2,2))) - 2TY_qc123(Z1 + Z2)pow(Z2,2))
821	t2 = ((Z1 + Z2)exp(Z1)pow(TY_qc223,2)pow(Z2,2) + (TY_q22 + TY_qc223)exp(2Z2)t1 + 2TY_qc223exp(Z2)(TY_qc123Z1pow(Z2,2) - (Z1 + Z2)exp(Z1)(TY_qa23 + (TY_qa23 + TY_qb23)Z2 + (TY_q22 + TY_qc223)*pow(Z2,2))))
822
823	TY_B25 = -6phiexp(-Z1 - 2Z2)pow(TY_q22,-2)pow(Z2,-3)t2*pow(Z1 + Z2,-1);
824
825
826	t1 = (TY_qa21TY_qc232 + TY_qa32TY_qc221(1 + Z2) + Z2(TY_qb32TY_qc221 + TY_qc232(TY_qa21 + TY_qb21 + 2TY_qc221Z2)))
827	t2 = (-(TY_qc132TY_qc221) - TY_qc121TY_qc232)
828	t3 = (TY_qb32TY_qc221(-1 + Z2)(Z1 + Z2) + TY_qb21TY_qc232(-1 + Z2)(Z1 + Z2) - Z2(TY_qc221TY_qc232Z1 + TY_qc121TY_qc232Z2 + TY_qc221(TY_q22 + TY_qc132 + TY_qc232)*Z2))
829	t4 = (exp(Z1)(2Z2t3 + TY_qa32TY_qc221(Z1 + Z2)(-2 + pow(Z2,2)) + TY_qa21TY_qc232(Z1 + Z2)(-2 + pow(Z2,2))) + 2(TY_qc132TY_qc221 + TY_qc121TY_qc232)(Z1 + Z2)pow(Z2,2))
830
831	TY_B32 = 6phiexp(-Z1 - 2Z2)pow(TY_q22,-2)pow(Z2,-3)(-2TY_qc221TY_qc232(Z1 + Z2)exp(Z1)pow(Z2,2) + 2exp(Z2)((Z1 + Z2)t1exp(Z1) + t2Z1pow(Z2,2)) + exp(2Z2)t4)pow(Z1 + Z2,-1);
832
833
834	t1 = (-((Z1 + Z2)(TY_qa23TY_qc232 + TY_qa32TY_qc223(1 + Z2) + Z2(TY_qb32TY_qc223 + TY_qc232(TY_qa23 + TY_qb23 + TY_q22Z2 + 2TY_qc223Z2)))exp(Z1)) + (TY_qc132TY_qc223 + TY_qc123TY_qc232)Z1*pow(Z2,2))
835	t2 = (TY_qc132(TY_q22 + TY_qc223) + TY_qc123TY_qc232)
836	t3 = (TY_qa32(TY_q22 + TY_qc223) + TY_qa23TY_qc232)
837	t4 = (TY_qa32(TY_q22 + TY_qc223) + TY_qa23TY_qc232 + (TY_qb32(TY_q22 + TY_qc223) + TY_qb23TY_qc232)*Z1)
838	t5 = (-2(TY_qb32(TY_q22 + TY_qc223) + TY_qb23TY_qc232) + ((TY_qa32 + 2TY_qb32)(TY_q22 + TY_qc223) + (-2TY_q22 + TY_qa23 + 2TY_qb23 - 2TY_qc223)TY_qc232)Z1)
839	t6 = (2t3Z1 + 2t4Z2 - t5pow(Z2,2) + ((2TY_q22 - TY_qa32 - 2TY_qb32 + 2TY_qc132)(TY_q22 + TY_qc223) + (2TY_q22 - TY_qa23 + 2(-TY_qb23 + TY_qc123 + TY_qc223))TY_qc232)*pow(Z2,3))
840
841	TY_B34 = -6phiexp(-Z1 - 2Z2)pow(TY_q22,-2)pow(Z2,-3)(2TY_qc223TY_qc232(Z1 + Z2)exp(Z1)pow(Z2,2) + 2exp(Z2)t1 + exp(2Z2)(-2t2(Z1 + Z2)pow(Z2,2) + exp(Z1)t6))pow(Z1 + Z2,-1);
842
843
844	// /*double norm_B = sqrt(pow(TY_B12, 2)+pow(TY_B14, 2)+pow(TY_B21, 2)+pow(TY_B22, 2)+pow(TY_B23, 2)+pow(TY_B24, 2)+pow(TY_B25, 2)+pow(TY_B32, 2)+pow(TY_B34, 2));
845	//
846	// TY_B12 /= norm_B;
847	// TY_B14 /= norm_B;
848	// TY_B21 /= norm_B;
849	// TY_B22 /= norm_B;
850	// TY_B23 /= norm_B;
851	// TY_B24 /= norm_B;
852	// TY_B25 /= norm_B;
853	// TY_B32 /= norm_B;
854	// TY_B34 /= norm_B; */
855
856	if( prnt )
857	printf "\rNonlinear equation 2 = \r"
858	// printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0., TY_B12, 0., TY_B14, 0.
859	// printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\r", TY_B21, TY_B22, TY_B23, TY_B24, TY_B25
860	// printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0., TY_B32, 0., TY_B34, 0.
861	printf "TY_B12 = %15.12g\r",TY_B12
862	printf "TY_B14 = %15.12g\r",TY_B14
863	printf "TY_B21 = %15.12g\r",TY_B21
864	printf "TY_B22 = %15.12g\r",TY_B22
865	printf "TY_B23 = %15.12g\r",TY_B23
866	printf "TY_B24 = %15.12g\r",TY_B24
867	printf "TY_B25 = %15.12g\r",TY_B25
868	printf "TY_B32 = %15.12g\r",TY_B32
869	printf "TY_B34 = %15.12g\r",TY_B34
870	endif
871
872	// /* decrease order of nonlinear equation 1 by means of equation 2 */
873	NVAR TY_F14 = root:yuk:TY_F14
874	NVAR TY_F16 = root:yuk:TY_F16
875	NVAR TY_F18 = root:yuk:TY_F18
876	NVAR TY_F23 = root:yuk:TY_F23
877	NVAR TY_F24 = root:yuk:TY_F24
878	NVAR TY_F25 = root:yuk:TY_F25
879	NVAR TY_F26 = root:yuk:TY_F26
880	NVAR TY_F27 = root:yuk:TY_F27
881	NVAR TY_F28 = root:yuk:TY_F28
882	NVAR TY_F29 = root:yuk:TY_F29
883	NVAR TY_F32 = root:yuk:TY_F32
884	NVAR TY_F33 = root:yuk:TY_F33
885	NVAR TY_F34 = root:yuk:TY_F34
886	NVAR TY_F35 = root:yuk:TY_F35
887	NVAR TY_F36 = root:yuk:TY_F36
888	NVAR TY_F37 = root:yuk:TY_F37
889	NVAR TY_F38 = root:yuk:TY_F38
890	NVAR TY_F39 = root:yuk:TY_F39
891	NVAR TY_F310 = root:yuk:TY_F310
892
893	TY_F14 = -(TY_A32TY_B12TY_B32) + TY_A52pow(TY_B12,2) + TY_A12pow(TY_B32,2)
894	TY_F16 = 2TY_A52TY_B12TY_B14 - TY_A32TY_B14TY_B32 - TY_A32TY_B12TY_B34 + 2TY_A12TY_B32TY_B34
895	TY_F18 = -(TY_A32TY_B14TY_B34) + TY_A52pow(TY_B14,2) + TY_A12pow(TY_B34,2)
896	TY_F23 = 2TY_A52TY_B12TY_B21 - TY_A41TY_B12TY_B32 - TY_A32TY_B21TY_B32 + TY_A21pow(TY_B32,2)
897	TY_F24 = 2TY_A52TY_B12TY_B22 - TY_A42TY_B12TY_B32 - TY_A32TY_B22TY_B32 + TY_A22pow(TY_B32,2)
898	TY_F25 = 2TY_A52TY_B14TY_B21 + 2TY_A52TY_B12TY_B23 - TY_A43TY_B12TY_B32 - TY_A41TY_B14TY_B32 - TY_A32TY_B23TY_B32 - TY_A41TY_B12TY_B34 - TY_A32TY_B21TY_B34 + 2TY_A21TY_B32TY_B34 + TY_A23pow(TY_B32,2)
899	TY_F26 = 2TY_A52TY_B14TY_B22 + 2TY_A52TY_B12TY_B24 - TY_A42TY_B14TY_B32 - TY_A32TY_B24TY_B32 - TY_A42TY_B12TY_B34 - TY_A32TY_B22TY_B34 + 2TY_A22TY_B32*TY_B34
900	TY_F27 = 2TY_A52TY_B14TY_B23 + 2TY_A52TY_B12TY_B25 - TY_A43TY_B14TY_B32 - TY_A32TY_B25TY_B32 - TY_A43TY_B12TY_B34 - TY_A41TY_B14TY_B34 - TY_A32TY_B23TY_B34 + 2TY_A23TY_B32TY_B34 + TY_A21pow(TY_B34,2)
901	TY_F28 = 2TY_A52TY_B14TY_B24 - TY_A42TY_B14TY_B34 - TY_A32TY_B24TY_B34 + TY_A22pow(TY_B34,2)
902	TY_F29 = 2TY_A52TY_B14TY_B25 - TY_A43TY_B14TY_B34 - TY_A32TY_B25TY_B34 + TY_A23pow(TY_B34,2)
903	TY_F32 = -(TY_A41TY_B21TY_B32) + TY_A52*pow(TY_B21,2)
904	TY_F33 = 2TY_A52TY_B21TY_B22 - TY_A42TY_B21TY_B32 - TY_A41TY_B22*TY_B32
905	TY_F34 = 2TY_A52TY_B21TY_B23 - TY_A43TY_B21TY_B32 - TY_A42TY_B22TY_B32 - TY_A41TY_B23TY_B32 - TY_A41TY_B21TY_B34 + TY_A52pow(TY_B22,2)
906	TY_F35 = 2TY_A52TY_B22TY_B23 + 2TY_A52TY_B21TY_B24 - TY_A43TY_B22TY_B32 - TY_A42TY_B23TY_B32 - TY_A41TY_B24TY_B32 - TY_A42TY_B21TY_B34 - TY_A41TY_B22TY_B34
907	TY_F36 = 2TY_A52TY_B22TY_B24 + 2TY_A52TY_B21TY_B25 - TY_A43TY_B23TY_B32 - TY_A42TY_B24TY_B32 - TY_A41TY_B25TY_B32 - TY_A43TY_B21TY_B34 - TY_A42TY_B22TY_B34 - TY_A41TY_B23TY_B34 + TY_A52*pow(TY_B23,2)
908	TY_F37 = 2TY_A52TY_B23TY_B24 + 2TY_A52TY_B22TY_B25 - TY_A43TY_B24TY_B32 - TY_A42TY_B25TY_B32 - TY_A43TY_B22TY_B34 - TY_A42TY_B23TY_B34 - TY_A41TY_B24TY_B34
909	TY_F38 = 2TY_A52TY_B23TY_B25 - TY_A43TY_B25TY_B32 - TY_A43TY_B23TY_B34 - TY_A42TY_B24TY_B34 - TY_A41TY_B25TY_B34 + TY_A52pow(TY_B24,2)
910	TY_F39 = 2TY_A52TY_B24TY_B25 - TY_A43TY_B24TY_B34 - TY_A42TY_B25*TY_B34
911	TY_F310 = -(TY_A43TY_B25TY_B34) + TY_A52*pow(TY_B25,2)
912
913	if( prnt )
914	printf "\rF = \r"
915	// printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0., 0., 0., TY_F14, 0., TY_F16, 0., TY_F18, 0., 0.
916	// printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0., 0., TY_F23, TY_F24, TY_F25, TY_F26, TY_F27, TY_F28, TY_F29, 0.
917	// printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0., TY_F32, TY_F33, TY_F34, TY_F35, TY_F36, TY_F37, TY_F38, TY_F39, TY_F310
918	printf "TY_F14 = %15.12g\r",TY_F14
919	printf "TY_F16 = %15.12g\r",TY_F16
920	printf "TY_F18 = %15.12g\r",TY_F18
921	printf "TY_F23 = %15.12g\r",TY_F23
922	printf "TY_F24 = %15.12g\r",TY_F24
923	printf "TY_F25 = %15.12g\r",TY_F25
924	printf "TY_F26 = %15.12g\r",TY_F26
925	printf "TY_F27 = %15.12g\r",TY_F27
926	printf "TY_F28 = %15.12g\r",TY_F28
927	printf "TY_F29 = %15.12g\r",TY_F29
928	printf "TY_F32 = %15.12g\r",TY_F32
929	printf "TY_F33 = %15.12g\r",TY_F33
930	printf "TY_F34 = %15.12g\r",TY_F34
931	printf "TY_F35 = %15.12g\r",TY_F35
932	printf "TY_F36 = %15.12g\r",TY_F36
933	printf "TY_F37 = %15.12g\r",TY_F37
934	printf "TY_F38 = %15.12g\r",TY_F38
935	printf "TY_F39 = %15.12g\r",TY_F39
936	printf "TY_F310 = %15.12g\r",TY_F310
937	endif
938
939	NVAR TY_G13 = root:yuk:TY_G13
940	NVAR TY_G14 = root:yuk:TY_G14
941	NVAR TY_G15 = root:yuk:TY_G15
942	NVAR TY_G16 = root:yuk:TY_G16
943	NVAR TY_G17 = root:yuk:TY_G17
944	NVAR TY_G18 = root:yuk:TY_G18
945	NVAR TY_G19 = root:yuk:TY_G19
946	NVAR TY_G110 = root:yuk:TY_G110
947	NVAR TY_G111 = root:yuk:TY_G111
948	NVAR TY_G112 = root:yuk:TY_G112
949	NVAR TY_G113 = root:yuk:TY_G113
950	NVAR TY_G22 = root:yuk:TY_G22
951	NVAR TY_G23 = root:yuk:TY_G23
952	NVAR TY_G24 = root:yuk:TY_G24
953	NVAR TY_G25 = root:yuk:TY_G25
954	NVAR TY_G26 = root:yuk:TY_G26
955	NVAR TY_G27 = root:yuk:TY_G27
956	NVAR TY_G28 = root:yuk:TY_G28
957	NVAR TY_G29 = root:yuk:TY_G29
958	NVAR TY_G210 = root:yuk:TY_G210
959	NVAR TY_G211 = root:yuk:TY_G211
960	NVAR TY_G212 = root:yuk:TY_G212
961	NVAR TY_G213 = root:yuk:TY_G213
962	NVAR TY_G214 = root:yuk:TY_G214
963
964
965	TY_G13 = -(TY_B12*TY_F32)
966	TY_G14 = -(TY_B12*TY_F33)
967	TY_G15 = TY_B32TY_F14 - TY_B14TY_F32 - TY_B12*TY_F34
968	TY_G16 = -(TY_B14TY_F33) - TY_B12TY_F35
969	TY_G17 = TY_B34TY_F14 + TY_B32TY_F16 - TY_B14TY_F34 - TY_B12TY_F36
970	TY_G18 = -(TY_B14TY_F35) - TY_B12TY_F37
971	TY_G19 = TY_B34TY_F16 + TY_B32TY_F18 - TY_B14TY_F36 - TY_B12TY_F38
972	TY_G110 = -(TY_B14TY_F37) - TY_B12TY_F39
973	TY_G111 = TY_B34TY_F18 - TY_B12TY_F310 - TY_B14*TY_F38
974	TY_G112 = -(TY_B14*TY_F39)
975	TY_G113 = -(TY_B14*TY_F310)
976	TY_G22 = -(TY_B21*TY_F32)
977	TY_G23 = -(TY_B22TY_F32) - TY_B21TY_F33
978	TY_G24 = TY_B32TY_F23 - TY_B23TY_F32 - TY_B22TY_F33 - TY_B21TY_F34
979	TY_G25 = TY_B32TY_F24 - TY_B24TY_F32 - TY_B23TY_F33 - TY_B22TY_F34 - TY_B21*TY_F35
980	TY_G26 = TY_B34TY_F23 + TY_B32TY_F25 - TY_B25TY_F32 - TY_B24TY_F33 - TY_B23TY_F34 - TY_B22TY_F35 - TY_B21*TY_F36
981	TY_G27 = TY_B34TY_F24 + TY_B32TY_F26 - TY_B25TY_F33 - TY_B24TY_F34 - TY_B23TY_F35 - TY_B22TY_F36 - TY_B21*TY_F37
982	TY_G28 = TY_B34TY_F25 + TY_B32TY_F27 - TY_B25TY_F34 - TY_B24TY_F35 - TY_B23TY_F36 - TY_B22TY_F37 - TY_B21*TY_F38
983	TY_G29 = TY_B34TY_F26 + TY_B32TY_F28 - TY_B25TY_F35 - TY_B24TY_F36 - TY_B23TY_F37 - TY_B22TY_F38 - TY_B21*TY_F39
984	TY_G210 = TY_B34TY_F27 + TY_B32TY_F29 - TY_B21TY_F310 - TY_B25TY_F36 - TY_B24TY_F37 - TY_B23TY_F38 - TY_B22*TY_F39
985	TY_G211 = TY_B34TY_F28 - TY_B22TY_F310 - TY_B25TY_F37 - TY_B24TY_F38 - TY_B23*TY_F39
986	TY_G212 = TY_B34TY_F29 - TY_B23TY_F310 - TY_B25TY_F38 - TY_B24TY_F39
987	TY_G213 = -(TY_B24TY_F310) - TY_B25TY_F39
988	TY_G214 = -(TY_B25*TY_F310)
989
990	if( prnt )
991	printf "\rG = \r"
992	// printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0., 0., TY_G13, TY_G14, TY_G15, TY_G16, TY_G17, TY_G18, TY_G19, TY_G110, TY_G111, TY_G112, TY_G113, 0.
993	// printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0., TY_G22, TY_G23, TY_G24, TY_G25, TY_G26, TY_G27, TY_G28, TY_G29, TY_G210, TY_G211, TY_G212, TY_G213, TY_G214
994	printf "TY_G13 = %15.12g\r",TY_G13
995	printf "TY_G14 = %15.12g\r",TY_G14
996	printf "TY_G15 = %15.12g\r",TY_G15
997	printf "TY_G16 = %15.12g\r",TY_G16
998	printf "TY_G17 = %15.12g\r",TY_G17
999	printf "TY_G18 = %15.12g\r",TY_G18
1000	printf "TY_G19 = %15.12g\r",TY_G19
1001	printf "TY_G110 = %15.12g\r",TY_G110
1002	printf "TY_G111 = %15.12g\r",TY_G111
1003	printf "TY_G112 = %15.12g\r",TY_G112
1004	printf "TY_G113 = %15.12g\r",TY_G113
1005	printf "TY_G22 = %15.12g\r",TY_G22
1006	printf "TY_G23 = %15.12g\r",TY_G23
1007	printf "TY_G24 = %15.12g\r",TY_G24
1008	printf "TY_G25 = %15.12g\r",TY_G25
1009	printf "TY_G26 = %15.12g\r",TY_G26
1010	printf "TY_G27 = %15.12g\r",TY_G27
1011	printf "TY_G28 = %15.12g\r",TY_G28
1012	printf "TY_G29 = %15.12g\r",TY_G29
1013	printf "TY_G210 = %15.12g\r",TY_G210
1014	printf "TY_G211 = %15.12g\r",TY_G211
1015	printf "TY_G212 = %15.12g\r",TY_G212
1016	printf "TY_G213 = %15.12g\r",TY_G213
1017	printf "TY_G214 = %15.12g\r",TY_G214
1018	endif
1019
1020	Make/O/D/N=23 TY_w
1021
1022	// coefficients for polynomial
1023	TY_w[0] = (-(TY_A21TY_B12) + TY_A12TY_B21)(TY_A52TY_B21 - TY_A41TY_B32)pow(TY_B21,2)*pow(TY_B32,3)
1024
1025	TY_w[1] = 2TY_B32TY_G13TY_G14 - TY_B24TY_G13TY_G22 - TY_B23TY_G14TY_G22 - TY_B22TY_G15TY_G22 - TY_B21TY_G16TY_G22 - TY_B23TY_G13TY_G23 - TY_B22TY_G14*TY_G23
1026	TY_w[1] += - TY_B21TY_G15TY_G23 + 2TY_B14TY_G22TY_G23 - TY_B22TY_G13TY_G24 - TY_B21TY_G14TY_G24 + 2TY_B12TY_G23TY_G24 - TY_B21TY_G13TY_G25 + 2TY_B12TY_G22*TY_G25
1027
1028	TY_w[2] = -(TY_B25TY_G13TY_G22) - TY_B24TY_G14TY_G22 - TY_B23TY_G15TY_G22 - TY_B22TY_G16TY_G22 - TY_B21TY_G17TY_G22 - TY_B24TY_G13TY_G23 - TY_B23TY_G14TY_G23 - TY_B22TY_G15TY_G23 - TY_B21TY_G16TY_G23
1029	TY_w[2] += -TY_B23TY_G13TY_G24 - TY_B22TY_G14TY_G24 - TY_B21TY_G15TY_G24 + 2TY_B14TY_G22TY_G24 - TY_B22TY_G13TY_G25 - TY_B21TY_G14TY_G25 + 2TY_B12TY_G23TY_G25 - TY_B21TY_G13TY_G26 + 2TY_B12TY_G22*TY_G26
1030	TY_w[2] += +TY_B34pow(TY_G13,2) + TY_B32(2TY_G13TY_G15 + pow(TY_G14,2)) + TY_B14pow(TY_G23,2) + TY_B12pow(TY_G24,2)
1031
1032	TY_w[3] = 2TY_B34TY_G13TY_G14 + 2TY_B32(TY_G14TY_G15 + TY_G13TY_G16) - TY_B25TY_G14TY_G22 - TY_B24TY_G15TY_G22 - TY_B23TY_G16TY_G22 - TY_B22TY_G17TY_G22 - TY_B21TY_G18TY_G22 - TY_B25TY_G13*TY_G23
1033	TY_w[3] += -TY_B24TY_G14TY_G23 - TY_B23TY_G15TY_G23 - TY_B22TY_G16TY_G23 - TY_B21TY_G17TY_G23 - TY_B24TY_G13TY_G24 - TY_B23TY_G14TY_G24 - TY_B22TY_G15TY_G24 - TY_B21TY_G16TY_G24 + 2TY_B14TY_G23*TY_G24
1034	TY_w[3] += -TY_B23TY_G13TY_G25 - TY_B22TY_G14TY_G25 - TY_B21TY_G15TY_G25 + 2TY_B14TY_G22TY_G25 + 2TY_B12TY_G24TY_G25 - TY_B22TY_G13TY_G26 - TY_B21TY_G14TY_G26 + 2TY_B12TY_G23TY_G26 - TY_B21TY_G13*TY_G27
1035	TY_w[3] += 2TY_B12TY_G22*TY_G27;
1036
1037	TY_w[4] = -(TY_B25TY_G15TY_G22) - TY_B24TY_G16TY_G22 - TY_B23TY_G17TY_G22 - TY_B22TY_G18TY_G22 - TY_B21TY_G19TY_G22 - TY_B25TY_G14TY_G23 - TY_B24TY_G15TY_G23 - TY_B23TY_G16TY_G23 - TY_B22TY_G17TY_G23
1038	TY_w[4] += -TY_B21TY_G18TY_G23 - TY_B25TY_G13TY_G24 - TY_B24TY_G14TY_G24 - TY_B23TY_G15TY_G24 - TY_B22TY_G16TY_G24 - TY_B21TY_G17TY_G24 - TY_B24TY_G13TY_G25 - TY_B23TY_G14TY_G25 - TY_B22TY_G15TY_G25
1039	TY_w[4] += -TY_B21TY_G16TY_G25 + 2TY_B14TY_G23TY_G25 - TY_B23TY_G13TY_G26 - TY_B22TY_G14TY_G26 - TY_B21TY_G15TY_G26 + 2TY_B14TY_G22TY_G26 + 2TY_B12TY_G24TY_G26 - TY_B22TY_G13TY_G27 - TY_B21TY_G14*TY_G27
1040	TY_w[4] += 2TY_B12TY_G23TY_G27 - TY_B21TY_G13TY_G28 + 2TY_B12TY_G22TY_G28 + TY_B34(2TY_G13TY_G15 + pow(TY_G14,2)) + TY_B32(2TY_G14TY_G16 + 2TY_G13TY_G17 + pow(TY_G15,2)) + TY_B14*pow(TY_G24,2)
1041	TY_w[4] += TY_B12*pow(TY_G25,2)
1042
1043	TY_w[5] = 2TY_B34(TY_G14TY_G15 + TY_G13TY_G16) + 2TY_B32(TY_G15TY_G16 + TY_G14TY_G17 + TY_G13TY_G18) - TY_B21TY_G110TY_G22 - TY_B25TY_G16TY_G22 - TY_B24TY_G17TY_G22 - TY_B23TY_G18TY_G22 - TY_B22TY_G19*TY_G22
1044	TY_w[5] += -TY_B25TY_G15TY_G23 - TY_B24TY_G16TY_G23 - TY_B23TY_G17TY_G23 - TY_B22TY_G18TY_G23 - TY_B21TY_G19TY_G23 - TY_B25TY_G14TY_G24 - TY_B24TY_G15TY_G24 - TY_B23TY_G16TY_G24 - TY_B22TY_G17TY_G24
1045	TY_w[5] += -TY_B21TY_G18TY_G24 - TY_B25TY_G13TY_G25 - TY_B24TY_G14TY_G25 - TY_B23TY_G15TY_G25 - TY_B22TY_G16TY_G25 - TY_B21TY_G17TY_G25 + 2TY_B14TY_G24TY_G25 - TY_B24TY_G13TY_G26 - TY_B23TY_G14*TY_G26
1046	TY_w[5] += -TY_B22TY_G15TY_G26 - TY_B21TY_G16TY_G26 + 2TY_B14TY_G23TY_G26 + 2TY_B12TY_G25TY_G26 - TY_B23TY_G13TY_G27 - TY_B22TY_G14TY_G27 - TY_B21TY_G15TY_G27 + 2TY_B14TY_G22TY_G27 + 2TY_B12TY_G24TY_G27
1047	TY_w[5] += -TY_B22TY_G13TY_G28 - TY_B21TY_G14TY_G28 + 2TY_B12TY_G23TY_G28 - TY_B21TY_G13TY_G29 + 2TY_B12TY_G22TY_G29
1048
1049	TY_w[6] = -(TY_B22TY_G110TY_G22) - TY_B21TY_G111TY_G22 - TY_B25TY_G17TY_G22 - TY_B24TY_G18TY_G22 - TY_B23TY_G19TY_G22 + TY_G210(-(TY_B21TY_G13) + 2TY_B12TY_G22) - TY_B21TY_G110TY_G23 - TY_B25TY_G16TY_G23
1050	TY_w[6] += -TY_B24TY_G17TY_G23 - TY_B23TY_G18TY_G23 - TY_B22TY_G19TY_G23 - TY_B25TY_G15TY_G24 - TY_B24TY_G16TY_G24 - TY_B23TY_G17TY_G24 - TY_B22TY_G18TY_G24 - TY_B21TY_G19TY_G24 - TY_B25TY_G14TY_G25
1051	TY_w[6] += -TY_B24TY_G15TY_G25 - TY_B23TY_G16TY_G25 - TY_B22TY_G17TY_G25 - TY_B21TY_G18TY_G25 - TY_B25TY_G13TY_G26 - TY_B24TY_G14TY_G26 - TY_B23TY_G15TY_G26 - TY_B22TY_G16TY_G26 - TY_B21TY_G17TY_G26
1052	TY_w[6] += 2TY_B14TY_G24TY_G26 - TY_B24TY_G13TY_G27 - TY_B23TY_G14TY_G27 - TY_B22TY_G15TY_G27 - TY_B21TY_G16TY_G27 + 2TY_B14TY_G23TY_G27 + 2TY_B12TY_G25TY_G27 - TY_B23TY_G13TY_G28 - TY_B22TY_G14*TY_G28
1053	TY_w[6] += -TY_B21TY_G15TY_G28 + 2TY_B14TY_G22TY_G28 + 2TY_B12TY_G24TY_G28 - TY_B22TY_G13TY_G29 - TY_B21TY_G14TY_G29 + 2TY_B12TY_G23TY_G29 + TY_B34(2TY_G14TY_G16 + 2TY_G13TY_G17 + pow(TY_G15,2))
1054	TY_w[6] += TY_B32(2(TY_G15TY_G17 + TY_G14TY_G18 + TY_G13TY_G19) + pow(TY_G16,2)) + TY_B14pow(TY_G25,2) + TY_B12*pow(TY_G26,2)
1055
1056	TY_w[7] = 2TY_B34(TY_G15TY_G16 + TY_G14TY_G17 + TY_G13TY_G18) + 2TY_B32(TY_G110TY_G13 + TY_G16TY_G17 + TY_G15TY_G18 + TY_G14TY_G19) - TY_B22TY_G13TY_G210 - TY_B21TY_G14TY_G210 - TY_B23TY_G110*TY_G22
1057	TY_w[7] += -TY_B22TY_G111TY_G22 - TY_B21TY_G112TY_G22 - TY_B25TY_G18TY_G22 - TY_B24TY_G19TY_G22 + TY_G211(-(TY_B21TY_G13) + 2TY_B12TY_G22) - TY_B22TY_G110TY_G23 - TY_B21TY_G111TY_G23 - TY_B25TY_G17TY_G23
1058	TY_w[7] += -TY_B24TY_G18TY_G23 - TY_B23TY_G19TY_G23 + 2TY_B12TY_G210TY_G23 - TY_B21TY_G110TY_G24 - TY_B25TY_G16TY_G24 - TY_B24TY_G17TY_G24 - TY_B23TY_G18TY_G24 - TY_B22TY_G19TY_G24 - TY_B25TY_G15*TY_G25
1059	TY_w[7] += -TY_B24TY_G16TY_G25 - TY_B23TY_G17TY_G25 - TY_B22TY_G18TY_G25 - TY_B21TY_G19TY_G25 - TY_B25TY_G14TY_G26 - TY_B24TY_G15TY_G26 - TY_B23TY_G16TY_G26 - TY_B22TY_G17TY_G26 - TY_B21TY_G18TY_G26
1060	TY_w[7] += 2TY_B14TY_G25TY_G26 - TY_B25TY_G13TY_G27 - TY_B24TY_G14TY_G27 - TY_B23TY_G15TY_G27 - TY_B22TY_G16TY_G27 - TY_B21TY_G17TY_G27 + 2TY_B14TY_G24TY_G27 + 2TY_B12TY_G26TY_G27 - TY_B24TY_G13*TY_G28
1061	TY_w[7] += -TY_B23TY_G14TY_G28 - TY_B22TY_G15TY_G28 - TY_B21TY_G16TY_G28 + 2TY_B14TY_G23TY_G28 + 2TY_B12TY_G25TY_G28 - TY_B23TY_G13TY_G29 - TY_B22TY_G14TY_G29 - TY_B21TY_G15TY_G29 + 2TY_B14TY_G22*TY_G29
1062	TY_w[7] += 2TY_B12TY_G24*TY_G29
1063
1064	TY_w[8] = -(TY_B23TY_G13TY_G210) - TY_B22TY_G14TY_G210 - TY_B21TY_G15TY_G210 - TY_B22TY_G13TY_G211 - TY_B21TY_G14TY_G211 - TY_B21TY_G13TY_G212 - TY_B24TY_G110TY_G22 - TY_B23TY_G111TY_G22 - TY_B22TY_G112TY_G22
1065	TY_w[8] += -TY_B21TY_G113TY_G22 - TY_B25TY_G19TY_G22 + 2TY_B14TY_G210TY_G22 + 2TY_B12TY_G212TY_G22 - TY_B23TY_G110TY_G23 - TY_B22TY_G111TY_G23 - TY_B21TY_G112TY_G23 - TY_B25TY_G18TY_G23 - TY_B24TY_G19TY_G23
1066	TY_w[8] += 2TY_B12TY_G211TY_G23 - TY_B22TY_G110TY_G24 - TY_B21TY_G111TY_G24 - TY_B25TY_G17TY_G24 - TY_B24TY_G18TY_G24 - TY_B23TY_G19TY_G24 + 2TY_B12TY_G210TY_G24 - TY_B21TY_G110TY_G25 - TY_B25TY_G16TY_G25
1067	TY_w[8] += -TY_B24TY_G17TY_G25 - TY_B23TY_G18TY_G25 - TY_B22TY_G19TY_G25 - TY_B25TY_G15TY_G26 - TY_B24TY_G16TY_G26 - TY_B23TY_G17TY_G26 - TY_B22TY_G18TY_G26 - TY_B21TY_G19TY_G26 - TY_B25TY_G14TY_G27
1068	TY_w[8] += -TY_B24TY_G15TY_G27 - TY_B23TY_G16TY_G27 - TY_B22TY_G17TY_G27 - TY_B21TY_G18TY_G27 + 2TY_B14TY_G25TY_G27 - TY_B25TY_G13TY_G28 - TY_B24TY_G14TY_G28 - TY_B23TY_G15TY_G28 - TY_B22TY_G16*TY_G28
1069	TY_w[8] += -TY_B21TY_G17TY_G28 + 2TY_B14TY_G24TY_G28 + 2TY_B12TY_G26TY_G28 - TY_B24TY_G13TY_G29 - TY_B23TY_G14TY_G29 - TY_B22TY_G15TY_G29 - TY_B21TY_G16TY_G29 + 2TY_B14TY_G23TY_G29 + 2TY_B12TY_G25TY_G29
1070	TY_w[8] += TY_B34(2(TY_G15TY_G17 + TY_G14TY_G18 + TY_G13TY_G19) + pow(TY_G16,2)) + TY_B32(2(TY_G111TY_G13 + TY_G110TY_G14 + TY_G16TY_G18 + TY_G15TY_G19) + pow(TY_G17,2)) + TY_B14pow(TY_G26,2)
1071	TY_w[8] += TY_B12*pow(TY_G27,2)
1072
1073	TY_w[9] = 2TY_B34(TY_G110TY_G13 + TY_G16TY_G17 + TY_G15TY_G18 + TY_G14TY_G19) + 2TY_B32(TY_G112TY_G13 + TY_G111TY_G14 + TY_G110TY_G15 + TY_G17TY_G18 + TY_G16TY_G19) - TY_B24TY_G13TY_G210 - TY_B23TY_G14*TY_G210
1074	TY_w[9] += -TY_B22TY_G15TY_G210 - TY_B21TY_G16TY_G210 - TY_B23TY_G13TY_G211 - TY_B22TY_G14TY_G211 - TY_B21TY_G15TY_G211 - TY_B22TY_G13TY_G212 - TY_B21TY_G14TY_G212 - TY_B25TY_G110TY_G22 - TY_B24TY_G111TY_G22
1075	TY_w[9] += -TY_B23TY_G112TY_G22 - TY_B22TY_G113TY_G22 + 2TY_B14TY_G211TY_G22 + TY_G213(-(TY_B21TY_G13) + 2TY_B12TY_G22) - TY_B24TY_G110TY_G23 - TY_B23TY_G111TY_G23 - TY_B22TY_G112*TY_G23
1076	TY_w[9] += -TY_B21TY_G113TY_G23 - TY_B25TY_G19TY_G23 + 2TY_B14TY_G210TY_G23 + 2TY_B12TY_G212TY_G23 - TY_B23TY_G110TY_G24 - TY_B22TY_G111TY_G24 - TY_B21TY_G112TY_G24 - TY_B25TY_G18TY_G24 - TY_B24TY_G19TY_G24
1077	TY_w[9] += 2TY_B12TY_G211TY_G24 - TY_B22TY_G110TY_G25 - TY_B21TY_G111TY_G25 - TY_B25TY_G17TY_G25 - TY_B24TY_G18TY_G25 - TY_B23TY_G19TY_G25 + 2TY_B12TY_G210TY_G25 - TY_B21TY_G110TY_G26 - TY_B25TY_G16TY_G26
1078	TY_w[9] += -TY_B24TY_G17TY_G26 - TY_B23TY_G18TY_G26 - TY_B22TY_G19TY_G26 - TY_B25TY_G15TY_G27 - TY_B24TY_G16TY_G27 - TY_B23TY_G17TY_G27 - TY_B22TY_G18TY_G27 - TY_B21TY_G19TY_G27 + 2TY_B14TY_G26*TY_G27
1079	TY_w[9] += -TY_B25TY_G14TY_G28 - TY_B24TY_G15TY_G28 - TY_B23TY_G16TY_G28 - TY_B22TY_G17TY_G28 - TY_B21TY_G18TY_G28 + 2TY_B14TY_G25TY_G28 + 2TY_B12TY_G27TY_G28 - TY_B25TY_G13TY_G29 - TY_B24TY_G14TY_G29
1080	TY_w[9] += -TY_B23TY_G15TY_G29 - TY_B22TY_G16TY_G29 - TY_B21TY_G17TY_G29 + 2TY_B14TY_G24TY_G29 + 2TY_B12TY_G26TY_G29
1081
1082	TY_w[10] = -(TY_B25TY_G13TY_G210) - TY_B24TY_G14TY_G210 - TY_B23TY_G15TY_G210 - TY_B22TY_G16TY_G210 - TY_B21TY_G17TY_G210 - TY_B24TY_G13TY_G211 - TY_B23TY_G14TY_G211 - TY_B22TY_G15TY_G211 - TY_B21TY_G16TY_G211
1083	TY_w[10] += -TY_B23TY_G13TY_G212 - TY_B22TY_G14TY_G212 - TY_B21TY_G15TY_G212 - TY_B22TY_G13TY_G213 - TY_B21TY_G14TY_G213 - TY_B21TY_G13TY_G214 - TY_B25TY_G111TY_G22 - TY_B24TY_G112TY_G22 - TY_B23TY_G113TY_G22
1084	TY_w[10] += 2TY_B14TY_G212TY_G22 + 2TY_B12TY_G214TY_G22 - TY_B25TY_G110TY_G23 - TY_B24TY_G111TY_G23 - TY_B23TY_G112TY_G23 - TY_B22TY_G113TY_G23 + 2TY_B14TY_G211TY_G23 + 2TY_B12TY_G213TY_G23
1085	TY_w[10] += -TY_B24TY_G110TY_G24 - TY_B23TY_G111TY_G24 - TY_B22TY_G112TY_G24 - TY_B21TY_G113TY_G24 - TY_B25TY_G19TY_G24 + 2TY_B14TY_G210TY_G24 + 2TY_B12TY_G212TY_G24 - TY_B23TY_G110TY_G25
1086	TY_w[10] += -TY_B22TY_G111TY_G25 - TY_B21TY_G112TY_G25 - TY_B25TY_G18TY_G25 - TY_B24TY_G19TY_G25 + 2TY_B12TY_G211TY_G25 - TY_B22TY_G110TY_G26 - TY_B21TY_G111TY_G26 - TY_B25TY_G17TY_G26 - TY_B24TY_G18*TY_G26
1087	TY_w[10] += -TY_B23TY_G19TY_G26 + 2TY_B12TY_G210TY_G26 - TY_B21TY_G110TY_G27 - TY_B25TY_G16TY_G27 - TY_B24TY_G17TY_G27 - TY_B23TY_G18TY_G27 - TY_B22TY_G19TY_G27 - TY_B25TY_G15TY_G28 - TY_B24TY_G16*TY_G28
1088	TY_w[10] += -TY_B23TY_G17TY_G28 - TY_B22TY_G18TY_G28 - TY_B21TY_G19TY_G28 + 2TY_B14TY_G26TY_G28 - TY_B25TY_G14TY_G29 - TY_B24TY_G15TY_G29 - TY_B23TY_G16TY_G29 - TY_B22TY_G17TY_G29 - TY_B21TY_G18*TY_G29
1089	TY_w[10] += 2TY_B14TY_G25TY_G29 + 2TY_B12TY_G27TY_G29 + TY_B34(2(TY_G111TY_G13 + TY_G110TY_G14 + TY_G16TY_G18 + TY_G15TY_G19) + pow(TY_G17,2))
1090	TY_w[10] += TY_B32(2(TY_G113TY_G13 + TY_G112TY_G14 + TY_G111TY_G15 + TY_G110TY_G16 + TY_G17TY_G19) + pow(TY_G18,2)) + TY_B14pow(TY_G27,2) + TY_B12*pow(TY_G28,2)
1091
1092	TY_w[11] = 2TY_B34(TY_G112TY_G13 + TY_G111TY_G14 + TY_G110TY_G15 + TY_G17TY_G18 + TY_G16TY_G19) + 2TY_B32(TY_G113TY_G14 + TY_G112TY_G15 + TY_G111TY_G16 + TY_G110TY_G17 + TY_G18TY_G19) - TY_B25TY_G14TY_G210
1093	TY_w[11] += -TY_B24TY_G15TY_G210 - TY_B23TY_G16TY_G210 - TY_B22TY_G17TY_G210 - TY_B21TY_G18TY_G210 - TY_B25TY_G13TY_G211 - TY_B24TY_G14TY_G211 - TY_B23TY_G15TY_G211 - TY_B22TY_G16TY_G211 - TY_B21TY_G17TY_G211
1094	TY_w[11] += -TY_B24TY_G13TY_G212 - TY_B23TY_G14TY_G212 - TY_B22TY_G15TY_G212 - TY_B21TY_G16TY_G212 - TY_B23TY_G13TY_G213 - TY_B22TY_G14TY_G213 - TY_B21TY_G15TY_G213 - TY_B25TY_G112TY_G22 - TY_B24TY_G113TY_G22
1095	TY_w[11] += 2TY_B14TY_G213TY_G22 - TY_B25TY_G111TY_G23 - TY_B24TY_G112TY_G23 - TY_B23TY_G113TY_G23 + 2TY_B14TY_G212TY_G23 - TY_G214(TY_B22TY_G13 + TY_B21TY_G14 - 2TY_B12TY_G23) - TY_B25TY_G110*TY_G24
1096	TY_w[11] += -TY_B24TY_G111TY_G24 - TY_B23TY_G112TY_G24 - TY_B22TY_G113TY_G24 + 2TY_B14TY_G211TY_G24 + 2TY_B12TY_G213TY_G24 - TY_B24TY_G110TY_G25 - TY_B23TY_G111TY_G25 - TY_B22TY_G112TY_G25
1097	TY_w[11] += -TY_B21TY_G113TY_G25 - TY_B25TY_G19TY_G25 + 2TY_B14TY_G210TY_G25 + 2TY_B12TY_G212TY_G25 - TY_B23TY_G110TY_G26 - TY_B22TY_G111TY_G26 - TY_B21TY_G112TY_G26 - TY_B25TY_G18TY_G26 - TY_B24TY_G19TY_G26
1098	TY_w[11] += 2TY_B12TY_G211TY_G26 - TY_B22TY_G110TY_G27 - TY_B21TY_G111TY_G27 - TY_B25TY_G17TY_G27 - TY_B24TY_G18TY_G27 - TY_B23TY_G19TY_G27 + 2TY_B12TY_G210TY_G27 - TY_B21TY_G110TY_G28 - TY_B25TY_G16TY_G28
1099	TY_w[11] += -TY_B24TY_G17TY_G28 - TY_B23TY_G18TY_G28 - TY_B22TY_G19TY_G28 + 2TY_B14TY_G27TY_G28 - TY_B25TY_G15TY_G29 - TY_B24TY_G16TY_G29 - TY_B23TY_G17TY_G29 - TY_B22TY_G18TY_G29 - TY_B21TY_G19*TY_G29
1100	TY_w[11] += 2TY_B14TY_G26TY_G29 + 2TY_B12TY_G28TY_G29
1101
1102	TY_w[12] = -(TY_B25TY_G15TY_G210) - TY_B24TY_G16TY_G210 - TY_B23TY_G17TY_G210 - TY_B22TY_G18TY_G210 - TY_B21TY_G19TY_G210 - TY_B25TY_G14TY_G211 - TY_B24TY_G15TY_G211 - TY_B23TY_G16TY_G211 - TY_B22TY_G17TY_G211
1103	TY_w[12] += -TY_B21TY_G18TY_G211 - TY_B25TY_G13TY_G212 - TY_B24TY_G14TY_G212 - TY_B23TY_G15TY_G212 - TY_B22TY_G16TY_G212 - TY_B21TY_G17TY_G212 - TY_B24TY_G13TY_G213 - TY_B23TY_G14TY_G213 - TY_B22TY_G15TY_G213
1104	TY_w[12] += -TY_B21TY_G16TY_G213 - TY_B25TY_G113TY_G22 - TY_B25TY_G112TY_G23 - TY_B24TY_G113TY_G23 + 2TY_B14TY_G213TY_G23 - TY_B25TY_G111TY_G24 - TY_B24TY_G112TY_G24 - TY_B23TY_G113*TY_G24
1105	TY_w[12] += 2TY_B14TY_G212TY_G24 - TY_G214(TY_B23TY_G13 + TY_B22TY_G14 + TY_B21TY_G15 - 2TY_B14TY_G22 - 2TY_B12TY_G24) - TY_B25TY_G110TY_G25 - TY_B24TY_G111TY_G25 - TY_B23TY_G112*TY_G25
1106	TY_w[12] += -TY_B22TY_G113TY_G25 + 2TY_B14TY_G211TY_G25 + 2TY_B12TY_G213TY_G25 - TY_B24TY_G110TY_G26 - TY_B23TY_G111TY_G26 - TY_B22TY_G112TY_G26 - TY_B21TY_G113TY_G26 - TY_B25TY_G19TY_G26
1107	TY_w[12] += 2TY_B14TY_G210TY_G26 + 2TY_B12TY_G212TY_G26 - TY_B23TY_G110TY_G27 - TY_B22TY_G111TY_G27 - TY_B21TY_G112TY_G27 - TY_B25TY_G18TY_G27 - TY_B24TY_G19TY_G27 + 2TY_B12TY_G211*TY_G27
1108	TY_w[12] += -TY_B22TY_G110TY_G28 - TY_B21TY_G111TY_G28 - TY_B25TY_G17TY_G28 - TY_B24TY_G18TY_G28 - TY_B23TY_G19TY_G28 + 2TY_B12TY_G210TY_G28 - TY_B21TY_G110TY_G29 - TY_B25TY_G16TY_G29 - TY_B24TY_G17*TY_G29
1109	TY_w[12] += -TY_B23TY_G18TY_G29 - TY_B22TY_G19TY_G29 + 2TY_B14TY_G27TY_G29 + TY_B34(2(TY_G113TY_G13 + TY_G112TY_G14 + TY_G111TY_G15 + TY_G110TY_G16 + TY_G17TY_G19) + pow(TY_G18,2))
1110	TY_w[12] += TY_B32(2(TY_G113TY_G15 + TY_G112TY_G16 + TY_G111TY_G17 + TY_G110TY_G18) + pow(TY_G19,2)) + TY_B14pow(TY_G28,2) + TY_B12pow(TY_G29,2)
1111
1112	TY_w[13] = 2TY_B32(TY_G113TY_G16 + TY_G112TY_G17 + TY_G111TY_G18 + TY_G110TY_G19) + 2TY_B34(TY_G113TY_G14 + TY_G112TY_G15 + TY_G111TY_G16 + TY_G110TY_G17 + TY_G18TY_G19) - TY_B21TY_G110*TY_G210
1113	TY_w[13] += -TY_B25TY_G16TY_G210 - TY_B24TY_G17TY_G210 - TY_B23TY_G18TY_G210 - TY_B22TY_G19TY_G210 - TY_B25TY_G15TY_G211 - TY_B24TY_G16TY_G211 - TY_B23TY_G17TY_G211 - TY_B22TY_G18TY_G211 - TY_B21TY_G19TY_G211
1114	TY_w[13] += -TY_B25TY_G14TY_G212 - TY_B24TY_G15TY_G212 - TY_B23TY_G16TY_G212 - TY_B22TY_G17TY_G212 - TY_B21TY_G18TY_G212 - TY_B25TY_G13TY_G213 - TY_B24TY_G14TY_G213 - TY_B23TY_G15TY_G213 - TY_B22TY_G16TY_G213
1115	TY_w[13] += -TY_B21TY_G17TY_G213 - TY_B25TY_G113TY_G23 - TY_B25TY_G112TY_G24 - TY_B24TY_G113TY_G24 + 2TY_B14TY_G213TY_G24 - TY_B25TY_G111TY_G25 - TY_B24TY_G112TY_G25 - TY_B23TY_G113*TY_G25
1116	TY_w[13] += 2TY_B14TY_G212TY_G25 - TY_G214(TY_B24TY_G13 + TY_B23TY_G14 + TY_B22TY_G15 + TY_B21TY_G16 - 2TY_B14TY_G23 - 2TY_B12TY_G25) - TY_B25TY_G110TY_G26 - TY_B24TY_G111TY_G26 - TY_B23TY_G112TY_G26
1117	TY_w[13] += -TY_B22TY_G113TY_G26 + 2TY_B14TY_G211TY_G26 + 2TY_B12TY_G213TY_G26 - TY_B24TY_G110TY_G27 - TY_B23TY_G111TY_G27 - TY_B22TY_G112TY_G27 - TY_B21TY_G113TY_G27 - TY_B25TY_G19TY_G27
1118	TY_w[13] += 2TY_B14TY_G210TY_G27 + 2TY_B12TY_G212TY_G27 - TY_B23TY_G110TY_G28 - TY_B22TY_G111TY_G28 - TY_B21TY_G112TY_G28 - TY_B25TY_G18TY_G28 - TY_B24TY_G19TY_G28 + 2TY_B12TY_G211*TY_G28
1119	TY_w[13] += -TY_B22TY_G110TY_G29 - TY_B21TY_G111TY_G29 - TY_B25TY_G17TY_G29 - TY_B24TY_G18TY_G29 - TY_B23TY_G19TY_G29 + 2TY_B12TY_G210TY_G29 + 2TY_B14TY_G28TY_G29
1120
1121	TY_w[14] = -(TY_B22TY_G110TY_G210) - TY_B21TY_G111TY_G210 - TY_B25TY_G17TY_G210 - TY_B24TY_G18TY_G210 - TY_B23TY_G19TY_G210 - TY_B21TY_G110TY_G211 - TY_B25TY_G16TY_G211 - TY_B24TY_G17TY_G211
1122	TY_w[14] += -TY_B23TY_G18TY_G211 - TY_B22TY_G19TY_G211 - TY_B25TY_G15TY_G212 - TY_B24TY_G16TY_G212 - TY_B23TY_G17TY_G212 - TY_B22TY_G18TY_G212 - TY_B21TY_G19TY_G212 - TY_B25TY_G14TY_G213 - TY_B24TY_G15TY_G213
1123	TY_w[14] += -TY_B23TY_G16TY_G213 - TY_B22TY_G17TY_G213 - TY_B21TY_G18TY_G213 - TY_B25TY_G113TY_G24 - TY_B25TY_G112TY_G25 - TY_B24TY_G113TY_G25 + 2TY_B14TY_G213TY_G25 - TY_B25TY_G111TY_G26 - TY_B24TY_G112*TY_G26
1124	TY_w[14] += -TY_B23TY_G113TY_G26 + 2TY_B14TY_G212TY_G26 - TY_G214(TY_B25TY_G13 + TY_B24TY_G14 + TY_B23TY_G15 + TY_B22TY_G16 + TY_B21TY_G17 - 2TY_B14TY_G24 - 2TY_B12TY_G26) - TY_B25TY_G110*TY_G27
1125	TY_w[14] += -TY_B24TY_G111TY_G27 - TY_B23TY_G112TY_G27 - TY_B22TY_G113TY_G27 + 2TY_B14TY_G211TY_G27 + 2TY_B12TY_G213TY_G27 - TY_B24TY_G110TY_G28 - TY_B23TY_G111TY_G28 - TY_B22TY_G112TY_G28
1126	TY_w[14] += -TY_B21TY_G113TY_G28 - TY_B25TY_G19TY_G28 + 2TY_B14TY_G210TY_G28 + 2TY_B12TY_G212TY_G28 - TY_B23TY_G110TY_G29 - TY_B22TY_G111TY_G29 - TY_B21TY_G112TY_G29 - TY_B25TY_G18TY_G29 - TY_B24TY_G19TY_G29
1127	TY_w[14] += 2TY_B12TY_G211TY_G29 + TY_B32(2(TY_G113TY_G17 + TY_G112TY_G18 + TY_G111TY_G19) + pow(TY_G110,2))
1128	TY_w[14] += TY_B34(2(TY_G113TY_G15 + TY_G112TY_G16 + TY_G111TY_G17 + TY_G110TY_G18) + pow(TY_G19,2)) + TY_B12pow(TY_G210,2) + TY_B14pow(TY_G29,2)
1129
1130	TY_w[15] = 2TY_B34(TY_G113TY_G16 + TY_G112TY_G17 + TY_G111TY_G18 + TY_G110TY_G19) + 2TY_B32(TY_G110TY_G111 + TY_G113TY_G18 + TY_G112TY_G19) - TY_B23TY_G110TY_G210 - TY_B22TY_G111*TY_G210
1131	TY_w[15] += -TY_B21TY_G112TY_G210 - TY_B25TY_G18TY_G210 - TY_B24TY_G19TY_G210 - TY_B22TY_G110TY_G211 - TY_B21TY_G111TY_G211 - TY_B25TY_G17TY_G211 - TY_B24TY_G18TY_G211 - TY_B23TY_G19TY_G211
1132	TY_w[15] += 2TY_B12TY_G210TY_G211 - TY_B21TY_G110TY_G212 - TY_B25TY_G16TY_G212 - TY_B24TY_G17TY_G212 - TY_B23TY_G18TY_G212 - TY_B22TY_G19TY_G212 - TY_B25TY_G15TY_G213 - TY_B24TY_G16*TY_G213
1133	TY_w[15] += -TY_B23TY_G17TY_G213 - TY_B22TY_G18TY_G213 - TY_B21TY_G19TY_G213 - TY_B25TY_G113TY_G25 - TY_B25TY_G112TY_G26 - TY_B24TY_G113TY_G26 + 2TY_B14TY_G213TY_G26 - TY_B25TY_G111TY_G27 - TY_B24TY_G112*TY_G27
1134	TY_w[15] += -TY_B23TY_G113TY_G27 + 2TY_B14TY_G212TY_G27 - TY_G214(TY_B25TY_G14 + TY_B24TY_G15 + TY_B23TY_G16 + TY_B22TY_G17 + TY_B21TY_G18 - 2TY_B14TY_G25 - 2TY_B12TY_G27) - TY_B25TY_G110*TY_G28
1135	TY_w[15] += -TY_B24TY_G111TY_G28 - TY_B23TY_G112TY_G28 - TY_B22TY_G113TY_G28 + 2TY_B14TY_G211TY_G28 + 2TY_B12TY_G213TY_G28 - TY_B24TY_G110TY_G29 - TY_B23TY_G111TY_G29 - TY_B22TY_G112TY_G29
1136	TY_w[15] += -TY_B21TY_G113TY_G29 - TY_B25TY_G19TY_G29 + 2TY_B14TY_G210TY_G29 + 2TY_B12TY_G212TY_G29
1137
1138	TY_w[16] = -(TY_B24TY_G110TY_G210) - TY_B23TY_G111TY_G210 - TY_B22TY_G112TY_G210 - TY_B21TY_G113TY_G210 - TY_B25TY_G19TY_G210 - TY_B23TY_G110TY_G211 - TY_B22TY_G111TY_G211 - TY_B21TY_G112TY_G211
1139	TY_w[16] += -TY_B25TY_G18TY_G211 - TY_B24TY_G19TY_G211 - TY_B22TY_G110TY_G212 - TY_B21TY_G111TY_G212 - TY_B25TY_G17TY_G212 - TY_B24TY_G18TY_G212 - TY_B23TY_G19TY_G212 + 2TY_B12TY_G210*TY_G212
1140	TY_w[16] += -TY_B21TY_G110TY_G213 - TY_B25TY_G16TY_G213 - TY_B24TY_G17TY_G213 - TY_B23TY_G18TY_G213 - TY_B22TY_G19TY_G213 - TY_B25TY_G113TY_G26 - TY_B25TY_G112TY_G27 - TY_B24TY_G113TY_G27
1141	TY_w[16] += 2TY_B14TY_G213TY_G27 - TY_B25TY_G111TY_G28 - TY_B24TY_G112TY_G28 - TY_B23TY_G113TY_G28 + 2TY_B14TY_G212TY_G28
1142	TY_w[16] += -TY_G214(TY_B25TY_G15 + TY_B24TY_G16 + TY_B23TY_G17 + TY_B22TY_G18 + TY_B21TY_G19 - 2TY_B14TY_G26 - 2TY_B12TY_G28) - TY_B25TY_G110TY_G29 - TY_B24TY_G111TY_G29 - TY_B23TY_G112TY_G29
1143	TY_w[16] += -TY_B22TY_G113TY_G29 + 2TY_B14TY_G211TY_G29 + 2TY_B12TY_G213TY_G29 + TY_B34(2(TY_G113TY_G17 + TY_G112TY_G18 + TY_G111*TY_G19) + pow(TY_G110,2))
1144	TY_w[16] += TY_B32(2TY_G110TY_G112 + 2TY_G113TY_G19 + pow(TY_G111,2)) + TY_B14pow(TY_G210,2) + TY_B12*pow(TY_G211,2)
1145
1146	TY_w[17] = 2TY_B32(TY_G111TY_G112 + TY_G110TY_G113) + 2TY_B34(TY_G110TY_G111 + TY_G113TY_G18 + TY_G112TY_G19) - TY_B25TY_G110TY_G210 - TY_B24TY_G111TY_G210 - TY_B23TY_G112TY_G210 - TY_B22TY_G113*TY_G210
1147	TY_w[17] += -TY_B24TY_G110TY_G211 - TY_B23TY_G111TY_G211 - TY_B22TY_G112TY_G211 - TY_B21TY_G113TY_G211 - TY_B25TY_G19TY_G211 + 2TY_B14TY_G210TY_G211 - TY_B23TY_G110TY_G212 - TY_B22TY_G111*TY_G212
1148	TY_w[17] += -TY_B21TY_G112TY_G212 - TY_B25TY_G18TY_G212 - TY_B24TY_G19TY_G212 + 2TY_B12TY_G211TY_G212 - TY_B22TY_G110TY_G213 - TY_B21TY_G111TY_G213 - TY_B25TY_G17TY_G213 - TY_B24TY_G18*TY_G213
1149	TY_w[17] += -TY_B23TY_G19TY_G213 + 2TY_B12TY_G210TY_G213 - TY_B25TY_G113TY_G27 - TY_B25TY_G112TY_G28 - TY_B24TY_G113TY_G28 + 2TY_B14TY_G213TY_G28 - TY_B25TY_G111TY_G29 - TY_B24TY_G112TY_G29
1150	TY_w[17] += -TY_B23TY_G113TY_G29 + 2TY_B14TY_G212TY_G29 - TY_G214(TY_B21TY_G110 + TY_B25TY_G16 + TY_B24TY_G17 + TY_B23TY_G18 + TY_B22TY_G19 - 2TY_B14TY_G27 - 2TY_B12*TY_G29)
1151
1152	TY_w[18] = -(TY_B25TY_G111TY_G210) - TY_B24TY_G112TY_G210 - TY_B23TY_G113TY_G210 - TY_B25TY_G110TY_G211 - TY_B24TY_G111TY_G211 - TY_B23TY_G112TY_G211 - TY_B22TY_G113TY_G211 - TY_B24TY_G110TY_G212
1153	TY_w[18] += -TY_B23TY_G111TY_G212 - TY_B22TY_G112TY_G212 - TY_B21TY_G113TY_G212 - TY_B25TY_G19TY_G212 + 2TY_B14TY_G210TY_G212 - TY_B23TY_G110TY_G213 - TY_B22TY_G111TY_G213 - TY_B21TY_G112*TY_G213
1154	TY_w[18] += -TY_B25TY_G18TY_G213 - TY_B24TY_G19TY_G213 + 2TY_B12TY_G211TY_G213 - TY_B25TY_G113*TY_G28
1155	TY_w[18] += -TY_G214(TY_B22TY_G110 + TY_B21TY_G111 + TY_B25TY_G17 + TY_B24TY_G18 + TY_B23TY_G19 - 2TY_B12TY_G210 - 2TY_B14TY_G28) - TY_B25TY_G112TY_G29 - TY_B24TY_G113TY_G29 + 2TY_B14TY_G213*TY_G29
1156	TY_w[18] += TY_B34(2TY_G110TY_G112 + 2TY_G113TY_G19 + pow(TY_G111,2)) + TY_B32(2TY_G111TY_G113 + pow(TY_G112,2)) + TY_B14pow(TY_G211,2) + TY_B12pow(TY_G212,2)
1157
1158	TY_w[19] = 2TY_B32TY_G112TY_G113 + 2TY_B34(TY_G111TY_G112 + TY_G110TY_G113) - TY_B25TY_G112TY_G210 - TY_B24TY_G113TY_G210 - TY_B25TY_G111TY_G211 - TY_B24TY_G112TY_G211 - TY_B23TY_G113*TY_G211
1159	TY_w[19] += -TY_B25TY_G110TY_G212 - TY_B24TY_G111TY_G212 - TY_B23TY_G112TY_G212 - TY_B22TY_G113TY_G212 + 2TY_B14TY_G211TY_G212 - TY_B24TY_G110TY_G213 - TY_B23TY_G111TY_G213 - TY_B22TY_G112*TY_G213
1160	TY_w[19] += -TY_B21TY_G113TY_G213 - TY_B25TY_G19TY_G213 + 2TY_B14TY_G210TY_G213 + 2TY_B12TY_G212TY_G213 - TY_B25TY_G113TY_G29
1161	TY_w[19] += -TY_G214(TY_B23TY_G110 + TY_B22TY_G111 + TY_B21TY_G112 + TY_B25TY_G18 + TY_B24TY_G19 - 2TY_B12TY_G211 - 2TY_B14TY_G29)
1162
1163	TY_w[20] = -(TY_B25TY_G113TY_G210) - TY_B25TY_G112TY_G211 - TY_B24TY_G113TY_G211 - TY_B25TY_G111TY_G212 - TY_B24TY_G112TY_G212 - TY_B23TY_G113TY_G212 - TY_B25TY_G110TY_G213 - TY_B24TY_G111TY_G213
1164	TY_w[20] += -TY_B23TY_G112TY_G213 - TY_B22TY_G113TY_G213 + 2TY_B14TY_G211TY_G213 - (TY_B24TY_G110 + TY_B23TY_G111 + TY_B22TY_G112 + TY_B21TY_G113 + TY_B25TY_G19 - 2TY_B14TY_G210 - 2TY_B12TY_G212)*TY_G214
1165	TY_w[20] += TY_B34(2TY_G111TY_G113 + pow(TY_G112,2)) + TY_B32pow(TY_G113,2) + TY_B14pow(TY_G212,2) + TY_B12pow(TY_G213,2)
1166
1167	TY_w[21] = TY_B25(TY_A23TY_B14(-3TY_A52TY_B24TY_B25 + (2TY_A43TY_B24 + TY_A42TY_B25)TY_B34) + TY_B25(TY_A22TY_B14(-(TY_A52TY_B25) + TY_A43TY_B34) + TY_A12(4TY_A52TY_B24TY_B25 - (3TY_A43TY_B24 + TY_A42TY_B25)TY_B34)))pow(TY_B34,3)
1168
1169	TY_w[22] = (-(TY_A23TY_B14) + TY_A12TY_B25)(TY_A52TY_B25 - TY_A43TY_B34)pow(TY_B25,2)*pow(TY_B34,3)
1170
1171	if( prnt )
1172	printf "\rCoefficients of polynomial\r"
1173	variable i
1174	for ( i = 0; i < 23; i+=1 )
1175	printf "w[%d] = %g\r", i, TY_w[i]
1176	endfor
1177	printf "\r"
1178	endif
1179
1180	end
1181
1182	Function TY_capQ( d2 )
1183	Variable d2
1184
1185	NVAR TY_B32 = root:yuk:TY_B32
1186	NVAR TY_B34 = root:yuk:TY_B34
1187
1188	return d2 * TY_B32 + pow( d2, 3 ) * TY_B34
1189	end
1190
1191	Function TY_V( d2 )
1192	variable d2
1193
1194	NVAR TY_G13 = root:yuk:TY_G13
1195	NVAR TY_G14 = root:yuk:TY_G14
1196	NVAR TY_G15 = root:yuk:TY_G15
1197	NVAR TY_G16 = root:yuk:TY_G16
1198	NVAR TY_G17 = root:yuk:TY_G17
1199	NVAR TY_G18 = root:yuk:TY_G18
1200	NVAR TY_G19 = root:yuk:TY_G19
1201	NVAR TY_G110 = root:yuk:TY_G110
1202	NVAR TY_G111 = root:yuk:TY_G111
1203	NVAR TY_G112 = root:yuk:TY_G112
1204	NVAR TY_G113 = root:yuk:TY_G113
1205
1206	return -( pow( d2, 2 ) * TY_G13 + pow( d2, 3 ) * TY_G14 + pow( d2, 4 ) * TY_G15 + pow( d2, 5 ) * TY_G16 + pow( d2, 6 ) * TY_G17 + pow( d2, 7 ) * TY_G18 + pow( d2, 8 ) * TY_G19 + pow( d2, 9 ) * TY_G110 + pow( d2, 10 ) * TY_G111 + pow( d2, 11 ) * TY_G112 + pow( d2, 12 ) * TY_G113 )
1207	end
1208
1209	Function TY_capW( d2 )
1210	Variable d2
1211
1212	variable tmp
1213
1214	NVAR TY_G22 = root:yuk:TY_G22
1215	NVAR TY_G23 = root:yuk:TY_G23
1216	NVAR TY_G24 = root:yuk:TY_G24
1217	NVAR TY_G25 = root:yuk:TY_G25
1218	NVAR TY_G26 = root:yuk:TY_G26
1219	NVAR TY_G27 = root:yuk:TY_G27
1220	NVAR TY_G28 = root:yuk:TY_G28
1221	NVAR TY_G29 = root:yuk:TY_G29
1222	NVAR TY_G210 = root:yuk:TY_G210
1223	NVAR TY_G211 = root:yuk:TY_G211
1224	NVAR TY_G212 = root:yuk:TY_G212
1225	NVAR TY_G213 = root:yuk:TY_G213
1226	NVAR TY_G214 = root:yuk:TY_G214
1227
1228
1229	tmp = d2 * TY_G22 + pow( d2, 2 ) * TY_G23 + pow( d2, 3 ) * TY_G24 + pow( d2, 4 ) * TY_G25 + pow( d2, 5 ) * TY_G26
1230	tmp += pow( d2, 6 ) * TY_G27 + pow( d2, 7 ) * TY_G28 + pow( d2, 8 ) * TY_G29 + pow( d2, 9 ) * TY_G210
1231	tmp += pow( d2, 10 ) * TY_G211 + pow( d2, 11 ) * TY_G212 + pow( d2, 12 ) * TY_G213 + pow( d2, 13 ) * TY_G214
1232
1233	return tmp
1234	end
1235
1236	Function TY_X( d2 )
1237	Variable d2
1238
1239	return TY_V( d2 ) / TY_capW( d2 )
1240	end
1241
1242	// solve the linear system depending on d1, d2 using Cramer's rule
1243	//
1244	// a,b,c1,c2 are passed by reference and returned
1245	//
1246	Function TY_SolveLinearEquations( d1, d2, a, b, c1, c2)
1247	Variable d1, d2, &a, &b, &c1, &c2
1248
1249
1250	NVAR TY_q22 = root:yuk:TY_q22
1251	NVAR TY_qa12 = root:yuk:TY_qa12
1252	NVAR TY_qa21 = root:yuk:TY_qa21
1253	NVAR TY_qa22 = root:yuk:TY_qa22
1254	NVAR TY_qa23 = root:yuk:TY_qa23
1255	NVAR TY_qa32 = root:yuk:TY_qa32
1256
1257	NVAR TY_qb12 = root:yuk:TY_qb12
1258	NVAR TY_qb21 = root:yuk:TY_qb21
1259	NVAR TY_qb22 = root:yuk:TY_qb22
1260	NVAR TY_qb23 = root:yuk:TY_qb23
1261	NVAR TY_qb32 = root:yuk:TY_qb32
1262
1263	NVAR TY_qc112 = root:yuk:TY_qc112
1264	NVAR TY_qc121 = root:yuk:TY_qc121
1265	NVAR TY_qc122 = root:yuk:TY_qc122
1266	NVAR TY_qc123 = root:yuk:TY_qc123
1267	NVAR TY_qc132 = root:yuk:TY_qc132
1268
1269	NVAR TY_qc212 = root:yuk:TY_qc212
1270	NVAR TY_qc221 = root:yuk:TY_qc221
1271	NVAR TY_qc222 = root:yuk:TY_qc222
1272	NVAR TY_qc223 = root:yuk:TY_qc223
1273	NVAR TY_qc232 = root:yuk:TY_qc232
1274
1275	Variable det = TY_q22 * d1 * d2
1276	Variable det_a = TY_qa12 * d2 + TY_qa21 * d1 + TY_qa22 * d1 * d2 + TY_qa23 * d1 * pow( d2, 2 ) + TY_qa32 * pow( d1, 2 ) * d2
1277	Variable det_b = TY_qb12 * d2 + TY_qb21 * d1 + TY_qb22 * d1 * d2 + TY_qb23 * d1 * pow( d2, 2 ) + TY_qb32 * pow( d1, 2 ) * d2
1278	Variable det_c1 = TY_qc112 * d2 + TY_qc121 * d1 + TY_qc122 * d1 * d2 + TY_qc123 * d1 * pow( d2, 2 ) + TY_qc132 * pow( d1, 2 ) * d2
1279	Variable det_c2 = TY_qc212 * d2 + TY_qc221 * d1 + TY_qc222 * d1 * d2 + TY_qc223 * d1 * pow( d2, 2 ) + TY_qc232 * pow( d1, 2 ) * d2
1280
1281	a = det_a / det
1282	b = det_b / det
1283	c1 = det_c1 / det
1284	c2 = det_c2 / det
1285	end
1286
1287	//Solve the system of linear and nonlinear equations for given Zi, Ki, phi which gives at
1288	// most 22 solutions for the parameters a,b,ci,di. From the set of solutions choose the
1289	// physical one and return it.
1290	//
1291	//
1292	// a,b,c1,c2,d1,d2 are passed by reference and returned
1293	//
1294	Function TY_SolveEquations( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2, prnt )
1295	Variable Z1, Z2, K1, K2, phi, &a, &b, &c1, &c2, &d1, &d2, prnt
1296
1297
1298	// reduce system to a polynomial from which all solution are extracted
1299	// by doing that a lot of global background variables are set
1300	TY_ReduceNonlinearSystem( Z1, Z2, K1, K2, phi, prnt )
1301
1302	// the two coupled non-linear eqautions were reduced to a
1303	// 22nd order polynomial, the roots are give all possible solutions
1304	// for d2, than d1 can be computed by the function X
1305
1306	Make/O/D/N=23 real_coefficient,imag_coefficient
1307	Make/O/D/N=22 real_root,imag_root
1308
1309	//integer degree of polynomial
1310	variable degree = 22
1311	Variable i
1312
1313	WAVE TY_w = TY_w
1314
1315
1316	////
1317	// now I need to replace this solution with FindRoots/P to get the polynomial roots
1318	////
1319
1320	// vector of real and imaginary coefficients in order of INCREASING powers
1321	for ( i = 0; i <= degree; i+=1 )
1322	// the global variablw TY_w was set by TY_ReduceNonlinearSystem
1323	real_coefficient[i] = TY_w[i]
1324	// imag_coefficient[i] = 0.;
1325	endfor
1326
1327	// zrhqr(real_coefficient, degree, NR_r, NR_i);
1328
1329	FindRoots/P=real_coefficient
1330
1331	WAVE/C W_polyRoots = W_polyRoots
1332
1333	for(i=0; i<degree; i+=1)
1334	real_root[i] = real(W_polyRoots[i])
1335	imag_root[i] = imag(W_polyRoots[i])
1336	endfor
1337
1338	//end - NR solution of polynomial
1339
1340
1341	// show the result if in debug mode
1342	Variable x, y
1343	if ( prnt )
1344	for ( i = 0; i < degree; i+=1 )
1345	x = real_root[i]
1346	y = imag_root[i]
1347	if ( chop( y ) == 0 )
1348	printf "root(%d) = %g\r", i+1, x
1349	else
1350	printf "root(%d) = %g + %g i\r", i+1, x, y
1351	endif
1352	endfor
1353	printf "\r"
1354	endif
1355
1356
1357
1358	// select real roots and those satisfying Q(x) != 0 and W(x) != 0
1359	// Paper: Cluster formation in two-Yukawa Fluids, J. Chem. Phys. 122, 2005
1360	// The right set of (a, b, c1, c2, d1, d2) should have the following properties:
1361	// (1) a > 0
1362	// (2) d1, d2 are real
1363	// (3) vi/Ki > 0 <=> g(Zi) > 0
1364	// (4) if there is still more than root, calculate g(r) for each root
1365	// and g(r) of the correct root should have the minimum average value
1366	// inside the hardcore
1367	Variable var_a, var_b, var_c1, var_c2, var_d1, var_d2
1368	Make/O/D/N=22 sol_a, sol_b, sol_c1, sol_c2, sol_d1, sol_d2
1369
1370	Variable j = 0
1371	for ( i = 0; i < degree; i+=1 )
1372
1373	x = real_root[i]
1374	y = imag_root[i]
1375
1376	if ( chop( y ) == 0 && TY_capW( x ) != 0 && TY_capQ( x ) != 0 )
1377
1378	var_d1 = TY_X( x )
1379	var_d2 = x
1380
1381	// solution of linear system for given d1, d2 to obtain a,b,ci,di
1382	// var_a, var_b, var_c1, var_c2 passed by reference
1383	TY_SolveLinearEquations( var_d1, var_d2, var_a, var_b, var_c1, var_c2 )
1384
1385	// select physical solutions, for details check paper: "Cluster formation in
1386	// two-Yukawa fluids", J. Chem. Phys. 122 (2005)
1387	if ( var_a > 0 && TY_g( Z1, phi, Z1, Z2, var_a, var_b, var_c1, var_c2, var_d1, var_d2 ) > 0 && TY_g( Z2, phi, Z1, Z2, var_a, var_b, var_c1, var_c2, var_d1, var_d2 ) > 0 )
1388	sol_a[j] = var_a
1389	sol_b[j] = var_b
1390	sol_c1[j] = var_c1
1391	sol_c2[j] = var_c2
1392	sol_d1[j] = var_d1
1393	sol_d2[j] = var_d2
1394
1395	if ( prnt )
1396	Variable eq1 = chop( TY_LinearEquation_1( Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] ) )
1397	Variable eq2 = chop( TY_LinearEquation_2( Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] ) )
1398	Variable eq3 = chop( TY_LinearEquation_3( Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] ) )
1399	Variable eq4 = chop( TY_LinearEquation_4( Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] ) )
1400	Variable eq5 = chop( TY_NonlinearEquation_1( Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] ) )
1401	Variable eq6 = chop( TY_NonlinearEquation_2( Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] ) )
1402
1403	printf "solution[%d] = (%g, %g, %g, %g, %g, %g), ( eq == 0 ) = (%g, %g, %g, %g, %g, %g)\r", j, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j], eq1 , eq2, eq3, eq4, eq5, eq6
1404	endif
1405
1406	j+=1
1407	endif //var_a >0...
1408	endif //chop
1409	endfor
1410	// number remaining roots
1411	Variable n_roots = j
1412
1413	// if there is still more than one root left, than choose the one with the minimum
1414	// average value inside the hardcore
1415	if ( n_roots > 1 )
1416
1417	/////
1418	// it seems like this section should all be replaced in bulk with internal FFT code, rather than slow integration
1419	//
1420	// -- also, be sure to handle r=0, or the sum will always be INF
1421	////
1422
1423	// the number of q values should be a power of 2
1424	// in order to speed up the FFT
1425	/// int n = 1 << 14;
1426	Variable n=16384 //2^14 points
1427
1428	// the maximum q value should be large enough
1429	// to enable a reasoble approximation of g(r)
1430	variable qmax = 16 * 10 * 2 * pi
1431	Variable q, dq = qmax / ( n - 1 )
1432
1433	// step size for g(r)
1434	variable dr
1435
1436	// allocate memory for pair correlation function g(r)
1437	// and structure factor S(q)
1438	Make/O/D/N=(n) sq,gr //gr will be redimensioned!!
1439
1440	// loop over all remaining roots
1441	Variable minVal = 1e50 //a really big number
1442	Variable selected_root = 10
1443	Variable sumVal = 0
1444
1445	for ( j = 0; j < n_roots; j+=1)
1446
1447	// calculate structure factor at different q values
1448	for ( i = 0; i < n; i+=1)
1449
1450	q = dq * i
1451	sq[i] = SqTwoYukawa( q, Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] )
1452
1453	if(i<10 && prnt)
1454	printf "after SqTwoYukawa: s(q) = %g\r",sq[i]
1455	endif
1456
1457	endfor
1458
1459	// calculate pair correlation function for given
1460	// structure factor, g(r) is computed at values
1461	// r(i) = i * dr
1462
1463	// Yuk_SqToGr( phi, dq, sq, dr, gr, n )
1464
1465
1466	Yuk_SqToGr_FFT( phi, dq, sq, dr, gr, n )
1467
1468	// determine sum inside the hardcore
1469	// 0 =< r < 1 of the pair-correlation function
1470	sumVal = 0
1471	for (i = 0; i < floor( 1. / dr ); i+=1 )
1472
1473	sumVal += abs( gr[i] )
1474
1475	if(i<10 && prnt)
1476	printf "g(r) in core = %g\r",abs(gr[i])
1477	endif
1478
1479	endfor
1480
1481	if ( sumVal < minVal )
1482	minVal = sumVal
1483	selected_root = j
1484	endif
1485
1486	if(prnt)
1487	printf "min = %g sum = %g\r",minVal,sumVal
1488	endif
1489
1490	endfor
1491
1492
1493	// physical solution was found
1494	a = sol_a [ selected_root ] //sol_a [ selected_root ];
1495	b = sol_b [ selected_root ]
1496	c1 = sol_c1[ selected_root ]
1497	c2 = sol_c2[ selected_root ]
1498	d1 = sol_d1[ selected_root ]
1499	d2 = sol_d2[ selected_root ]
1500
1501	return 1
1502
1503	else
1504	if ( n_roots == 1 )
1505
1506	a = sol_a [0]
1507	b = sol_b [0]
1508	c1 = sol_c1[0]
1509	c2 = sol_c2[0]
1510	d1 = sol_d1[0]
1511	d2 = sol_d2[0]
1512
1513	return 1
1514	else
1515	// no solution was found
1516	return 0
1517	endif
1518	endif
1519
1520	end
1521
1522
1523
1524
1525	//
1526	Function Yuk_SqToGr_FFT( phi, dq, sq, dr, gr, n )
1527	Variable phi, dq
1528	WAVE sq
1529	Variable &dr
1530	WAVE gr
1531	Variable n
1532
1533	Variable npts,ii,rval,jj,qval,spread=1
1534	Variable alpha
1535
1536
1537	WaveStats/Q sq
1538	npts = V_npnts
1539
1540	dr = 2pi/(nptsdq)
1541
1542	Make/O/D/N=(npts) temp
1543
1544	temp = p*(sq[p] - 1)
1545	alpha = npts * pow( dq, 3 ) / ( 24 * pi * pi * phi )
1546
1547	FFT/OUT=1/DEST=W_FFT temp
1548
1549
1550	WAVE/C W_FFT = W_FFT
1551
1552	Redimension/N=(numpnts(W_FFT)) gr
1553
1554	gr = 1 + alpha/p*imag(W_FFT)
1555
1556	gr[0] = 0
1557
1558	// Killwaves/Z temp
1559
1560	return(0)
1561
1562	End
1563
1564	////////////////////////end converted procedures //////////////////////////////////

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format