Context Navigation

Two_Yukawa_v40.ipf @ 773

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

Changed the name of the Yukawa procedure to not end in _sq so that its help file could be found.

Updated the ModelDocs? to include the remaining 2010 model functions - Raspberry and Yukawa.

File size: 89.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	// -- 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
54	Proc 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","charge, 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")
75	End
76
77	//AAO version
78	Function 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)
91	End
92
93
94	// no igor code, return 0
95	//
96	Function fOneYukawa(w,x) : FitFunc
97	Wave w
98	Variable x
99
100	return (0)
101	End
102
103	//////////////////////////////////////////////////////////////
104	Proc 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","charge, Z1","scale, K2","charge, 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
127	End
128
129
130	//AAO version
131	Function 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)
154	End
155
156	Proc 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)
162	End
163
164	//
165	Function fTwoYukawa(cw,xw,yw) : FitFunc
166	Wave cw,xw,yw
167
168	Variable Z1, Z2, K1, K2, phi,radius
169	phi = cw[0]
170	radius = cw[1]
171	K1 = cw[2]
172	Z1 = cw[3]
173	K2 = cw[4]
174	Z2 = cw[5]
175
176	Variable a,b,c1,c2,d1,d2
177
178	Variable ok,check,prnt
179	prnt = 0 //print out intermediates
180
181	ok = TY_SolveEquations( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2, prnt ) // a,b,c1,c2,d1,d2 are returned
182	if(ok)
183	check = TY_CheckSolution( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
184	if(prnt)
185	printf "solution = (%g, %g, %g, %g, %g, %g) check = %d\r", a, b, c1, c2, d1, d2, check
186	endif
187
188	// if(check)
189	if(ok) //if(ok) simply takes the best solution, not necessarily one that passes TY_CheckSolution
190	yw = SqTwoYukawa(xwradius2, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2)
191	// printf("%g %g\n",q,sq)
192	endif
193	endif
194
195	return (0)
196	End
197
198
199	Macro Plot_2YukawaPotential()
200	fPlot_2YukawaPotential()
201	End
202
203	Function fPlot_2YukawaPotential()
204
205	Variable k1,z1,k2,z2,radius
206	Variable ii=0,num=500,rmax=10,rval
207
208	if(exists("root:coef_2yuk") == 0)
209	Abort "You must plot the 2-Yukawa model before plotting the potential"
210	else
211	WAVE coef_2yuk = root:coef_2yuk
212	endif
213
214	radius = coef_2yuk[1]
215	K1 = coef_2yuk[2]
216	Z1 = coef_2yuk[3]
217	K2 = coef_2yuk[4]
218	Z2 = coef_2yuk[5]
219
220	Make/O/D/N=(num) TwoYukawa_Potential,TwoYukawa_Potential_r
221	TwoYukawa_Potential_r = x/num*rmax
222
223	do
224	rval = TwoYukawa_Potential_r[ii]
225	if(rval <= 1)
226	TwoYukawa_Potential[ii] = inf
227	else
228	TwoYukawa_Potential[ii] = -1K1(exp(-1Z1(rval-1))/rval) - K2exp(-1Z2*(rval-1))/rval
229	endif
230
231	ii+=1
232	while(ii<num)
233
234
235	// if graph is not open, draw a graph
236	DoWindow YukawaPotential
237	if(V_flag == 0)
238	Display/N=YukawaPotential TwoYukawa_Potential vs TwoYukawa_Potential_r
239	ModifyGraph marker=29,msize=2,mode=4,grid=1,mirror=2
240	Label bottom "r/Diameter"
241	Label left "V(r)/kT"
242	endif
243
244	return(0)
245	End
246
247
248
249	///////////////////// converted procedures from c-code ////////////////////////////
250
251	/// 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
252
253	// function TY_W change to TY_capW, since there is a wave named TY_w
254
255
256
257
258
259	Static Function chop(x)
260	Variable x
261
262	if ( abs(x) < 1e-6 )
263	return 0
264	else
265	return x
266	endif
267
268	end
269
270	Static Function pow(a,b)
271	Variable a,b
272
273	return (a^b)
274	end
275
276	///*
277	// ==================================================================================================
278	//
279	// The two-yukawa structure factor is uniquley determined by 6 parameters a, b, c1, c2, d1, d2,
280	// which are the solution of a system of 6 equations ( 4 linear, 2 nonlinear ). The solution can
281	// constructed by the roots of a polynomial of 22nd degree. For more details see attached
282	// Mathematica notebook, where a derivation is given
283	//
284	// ==================================================================================================
285	// */
286
287	// these all may need to be declared as global variables !!
288	//
289	// - they are defined in a global scope in the c-code!
290	//
291	// - change the data folder
292	Function declare2YGlobals()
293
294	NewDataFolder/O/S root:yuk
295
296	Variable/G TY_q22
297	Variable/G TY_qa12, TY_qa21, TY_qa22, TY_qa23, TY_qa32
298	Variable/G TY_qb12, TY_qb21, TY_qb22, TY_qb23, TY_qb32
299	Variable/G TY_qc112, TY_qc121, TY_qc122, TY_qc123, TY_qc132
300	Variable/G TY_qc212, TY_qc221, TY_qc222, TY_qc223, TY_qc232
301	Variable/G TY_A12, TY_A21, TY_A22, TY_A23, TY_A32, TY_A41, TY_A42, TY_A43, TY_A52
302	Variable/G TY_B12, TY_B14, TY_B21, TY_B22, TY_B23, TY_B24, TY_B25, TY_B32, TY_B34
303	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
304	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
305
306	SetDataFolder root:
307	//this is an array, already global TY_w[23];
308
309	End
310
311
312	Function TY_sigma( s, Z1, Z2, a, b, c1, c2, d1, d2 )
313	Variable s, Z1, Z2, a, b, c1, c2, d1, d2
314
315	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 )
316	end
317
318	Function TY_tau( s, Z1, Z2, a, b, c1, c2 )
319	Variable s, Z1, Z2, a, b, c1, c2
320
321	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 ) )
322	end
323
324	Function TY_q( s, Z1, Z2, a, b, c1, c2, d1, d2 )
325	Variable s, Z1, Z2, a, b, c1, c2, d1, d2
326	return TY_sigma(s, Z1, Z2, a, b, c1, c2, d1, d2) - exp( -s ) * TY_tau(s, Z1, Z2, a,b, c1, c2)
327	end
328
329	Function TY_g( s, phi, Z1, Z2, a, b, c1, c2, d1, d2 )
330	Variable s, phi, Z1, Z2, a, b, c1, c2, d1, d2
331	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 ) )
332	end
333
334	///*
335	// ==================================================================================================
336	//
337	// Structure factor for the potential
338	//
339	// V(r) = -kB * T * ( K1 * exp[ -Z1 * (r - 1)] / r + K2 * exp[ -Z2 * (r - 1)] / r ) for r > 1
340	// V(r) = inf for r <= 1
341	//
342	// The structure factor is parametrized by (a, b, c1, c2, d1, d2)
343	// which depend on (K1, K2, Z1, Z2, phi).
344	//
345	// ==================================================================================================
346	// */
347
348	Function TY_hq( q, Z, K, v )
349	Variable q, Z, K, v
350
351	if ( q == 0)
352	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.
353	else
354
355	variable t1, t2, t3, t4
356
357	t1 = ( 1. - v / ( 2. * K * Z * exp( Z ) ) ) * ( ( 1. - cos( q ) ) / ( qq ) - 1. / ( ZZ + q*q ) )
358	t2 = ( vv ( q * cos( q ) - Z * sin( q ) ) ) / ( 4. * K * ZZ q * ( ZZ + qq ) )
359	t3 = ( q * cos( q ) + Z * sin( q ) ) / ( q * ( ZZ + qq ) )
360	t4 = v / ( Z * exp( Z ) ) - vv / ( 4. K * ZZ exp( 2. * Z ) ) - K
361
362	return v / Z * t1 - t2 + t3 * t4
363	endif
364	end
365
366
367	Function TY_pc( q, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
368	Variable q, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
369
370	variable v1 = 24. * phi * K1 * exp( Z1 ) * TY_g( Z1, phi, Z1, Z2, a, b, c1, c2, d1, d2 )
371	variable v2 = 24. * phi * K2 * exp( Z2 ) * TY_g( Z2, phi, Z1, Z2, a, b, c1, c2, d1, d2 )
372
373	variable a0 = a * a
374	variable b0 = -12. * phi ( pow( a + b,2 ) / 2. + a ( c1 * exp( -Z1 ) + c2 * exp( -Z2 ) ) )
375
376	variable t1, t2, t3
377
378	if ( q == 0 )
379	t1 = a0 / 3.
380	t2 = b0 / 4.
381	t3 = a0 * phi / 12.
382	else
383	t1 = a0 * ( sin( q ) - q * cos( q ) ) / pow( q, 3 )
384	t2 = b0 * ( 2. * q * sin( q ) - ( q * q - 2. ) * cos( q ) - 2. ) / pow( q, 4 )
385	t3 = a0 * phi * ( ( qq - 6. ) 4. * q * sin( q ) - ( pow( q, 4 ) - 12. * qq + 24.) cos( q ) + 24. ) / ( 2. * pow( q, 6 ) )
386	endif
387
388	variable t4 = TY_hq( q, Z1, K1, v1 ) + TY_hq( q, Z2, K2, v2 )
389
390	return -24. * phi * ( t1 + t2 + t3 + t4 )
391	end
392
393	Function SqTwoYukawa( q, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
394	variable q, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
395
396	if ( Z1 == Z2 )
397	// one-yukawa potential
398	return 0
399	else
400	// two-yukawa potential
401	return 1. / ( 1. - TY_pc( q, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
402	endif
403	end
404
405	///*
406	//==================================================================================================
407	//
408	// Non-linear eqaution system that determines the parameter for structure factor
409	//
410	//==================================================================================================
411	//*/
412
413	Function TY_LinearEquation_1( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
414	Variable Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
415
416	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 ) ) ) )
417	end
418
419	Function TY_LinearEquation_2( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
420	Variable Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
421
422	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 ) ) )
423	end
424
425	Function TY_LinearEquation_3( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
426	Variable Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
427
428	return K1 * exp( Z1 ) - d1 * Z1 * ( 1. - 12. * phi * TY_q( Z1, Z1, Z2, a, b, c1, c2, d1, d2 ) )
429	end
430
431	Function TY_LinearEquation_4( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
432	Variable Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
433
434	return K2 * exp( Z2 ) - d2 * Z2 * ( 1. - 12. * phi * TY_q( Z2, Z1, Z2, a, b, c1, c2, d1, d2 ) )
435	end
436
437	Function TY_NonlinearEquation_1( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
438	Variable Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
439
440	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 ) )
441	end
442
443	Function TY_NonlinearEquation_2( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
444	Variable Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
445
446	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 ) )
447	end
448
449	// Check the computed solutions satisfy the system of equations
450	Function TY_CheckSolution( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 )
451	variable Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2
452
453	variable eq_1 = chop( TY_LinearEquation_1( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
454	variable eq_2 = chop( TY_LinearEquation_2( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
455	variable eq_3 = chop( TY_LinearEquation_3( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
456	variable eq_4 = chop( TY_LinearEquation_4( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
457	variable eq_5 = chop( TY_NonlinearEquation_1( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
458	variable eq_6 = chop( TY_NonlinearEquation_2( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) )
459
460	// printf("Check of solution = %g %g %g %g %g %g\r",eq_1,eq_2,eq_3,eq_4,eq_5,eq_6);
461	// check if all equation are zero
462	return ( eq_1 == 0 && eq_2 == 0 && eq_3 == 0 && eq_4 == 0 && eq_5 == 0 && eq_6 == 0 )
463	end
464
465	Function TY_ReduceNonlinearSystem( Z1, Z2, K1, K2, phi, prnt )
466	Variable Z1, Z2, K1, K2, phi, prnt
467
468
469	// /* solution of the 4 linear equations depending on d1 and d2, the solution is polynomial
470	// in d1, d2. We represend the solution as determiants obtained by Cramer's rule
471	// which can be expressed by their coefficient matrices
472	// */
473
474	Variable m11 = (3.*phi)/2.
475	Variable m13 = 6.phiexp(-Z1)(2. + Z1(2. + Z1) - 2.exp(Z1))pow(Z1,-2)
476	Variable m14 = 6.phiexp(-Z2)(2. + Z2(2. + Z2) - 2.exp(Z2))pow(Z2,-2)
477	Variable m23 = -12.phiexp(-Z1)(-1. - Z1 + exp(Z1))pow(Z1,-1)
478	Variable m24 = -12.phiexp(-Z2)(-1. - Z2 + exp(Z2))pow(Z2,-1)
479	Variable m31 = -6.phiexp(-Z1)pow(Z1,-2)(2.(1 + Z1) + exp(Z1)(-2. + pow(Z1,2)))
480	Variable m32 = -12.phi(-1. + Z1 + exp(-Z1))*pow(Z1,-1)
481	Variable m33 = 6.phiexp(-2.Z1)pow(-1. + exp(Z1),2)
482	Variable m34 = 12.phiexp(-Z1 - Z2)(Z2 - (Z1 + Z2)exp(Z1) + Z1exp(Z1 + Z2))pow(Z1 + Z2,-1)
483	Variable m41 = -6.phiexp(-Z2)pow(Z2,-2)(2.(1. + Z2) + exp(Z2)(-2. + pow(Z2,2)))
484	Variable m42 = -12.phi(-1. + Z2 + exp(-Z2))*pow(Z2,-1)
485	Variable m43 = 12.phiexp(-Z1 - Z2)(Z1 - (Z1 + Z2 - Z2exp(Z1))exp(Z2))pow(Z1 + Z2,-1)
486	Variable m44 = 6.phiexp(-2Z2)pow(-1. + exp(Z2),2)
487
488	// /* determinant of the linear system expressed as coefficient matrix in d1, d2 */
489
490	NVAR TY_q22 = root:yuk:TY_q22
491
492	TY_q22 = m14(-(m33m42) + m23(m32m41 - m31m42) + m32m43 + (4.m11(-3.m33m41 + 2.m33m42 + 3.m31m43 - 2.m32m43))/3.)
493	TY_q22 += m13(m34m42 + m24(-(m32m41) + m31m42) - m32m44 + (4.m11(3.m34m41 - 2.m34m42 - 3.m31m44 + 2.m32m44))/3.)
494	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.
495
496	if( prnt )
497	printf "\rDet = \r"
498	// printf "%f\t%f\r%f\t%f\r", 0., 0., 0., TY_q22
499	printf "TY_q22 = %15.12g\r",TY_q22
500	endif
501
502	// /* Matrix representation of the determinant of the of the system where row refering to
503	// the variable a is replaced by solution vector */
504
505	NVAR TY_qa12 = root:yuk:TY_qa12
506	NVAR TY_qa21 = root:yuk:TY_qa21
507	NVAR TY_qa22 = root:yuk:TY_qa22
508	NVAR TY_qa23 = root:yuk:TY_qa23
509	NVAR TY_qa32 = root:yuk:TY_qa32
510
511	Variable t1,t2,t3,t4,t5,t6,t7,t8,t9,t10
512	Variable t11,t12,t13,t14,t15,t16,t17,t18,t19,t20 //simply to keep the line length small enough
513
514	TY_qa12 = (K1(3.m14(m23m42 - 4.m11m43) - 3.m13(m24m42 - 4.m11m44) + (3. + 4.m11)(m24m43 - m23m44))exp(Z1))/3.
515
516	TY_qa21 = -(K2(3.m14(m23m32 - 4.m11m33) - 3.m13(m24m32 - 4.m11m34) + (3. + 4.m11)(m24m33 - m23m34))exp(Z2))/3.
517
518	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.
519
520
521	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)))
522	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)))
523	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))))
524
525	TY_qa23 = 2.phipow(Z2,-2)(m24t1 + m23t2 + 2.t3)*pow(Z1 + Z2,-1)
526
527
528
529	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)))
530	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)))
531	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))))
532
533	TY_qa32 = 2.phipow(Z1,-2)(m24t1 + m23t2 + 2.t3)*pow(Z1 + Z2,-1)
534
535	if( prnt )
536	printf "\rDet_a = \r"
537	// 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.
538	printf "TY_qa12 = %15.12g\r",TY_qa12
539	printf "TY_qa21 = %15.12g\r",TY_qa21
540	printf "TY_qa22 = %15.12g\r",TY_qa22
541	printf "TY_qa23 = %15.12g\r",TY_qa23
542	printf "TY_qa32 = %15.12g\r",TY_qa32
543	endif
544
545	// /* Matrix representation of the determinant of the of the system where row refering to
546	// the variable b is replaced by solution vector */
547
548	NVAR TY_qb12 = root:yuk:TY_qb12
549	NVAR TY_qb21 = root:yuk:TY_qb21
550	NVAR TY_qb22 = root:yuk:TY_qb22
551	NVAR TY_qb23 = root:yuk:TY_qb23
552	NVAR TY_qb32 = root:yuk:TY_qb32
553
554	TY_qb12 = (K1(-3.m11m24m43 + m14(-3.m23m41 + (-3. + 8.m11)m43) + 3.m11m23m44 + m13(3.m24m41 + 3.m44 - 8.m11m44))*exp(Z1))/3.
555
556	TY_qb21 = (K2(-3.m13m24m31 + 3.m11m24m33 + m14(3.m23m31 + (3. - 8.m11)m33) - 3.m13m34 + 8.m11m13m34 - 3.m11m23m34)*exp(Z2))/3.
557
558	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))
559
560	t1 = (-(m14m33m41) + m13m34m41 + m14m31m43 - m11m34m43 - m13m31m44 + m11m33m44)
561	t2 = (-3.m11m24m43 + m14(-3.m23m41 + (-3. + 8.m11)m43) + 3.m11m23m44 + m13(3.m24m41 + 3.m44 - 8.m11*m44))
562	t3 = (3.m24(m33m41 - m31m43) + m23(-3.m34m41 + 3.m31m44) + (-3. + 8.m11)(m34m43 - m33*m44))
563
564	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))
565
566
567	t1 = (-(m34(m23m41 + m43)) + m24(m33m41 - m31m43) + (m23m31 + m33)*m44)
568	t2 = (-(m14m33m41) + m13m34m41 + m14m31m43 - m13m31m44)
569	t3 = (m14(2.m23m31 + 2.m33 - m23m41 - m43) + m13(-2.m34 + m24(-2.*m31 + m41) + m44))
570	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))
571	t5 = (2.m34m43(8. - 3.Z1) + 2.m33m44(-8. + 3.Z1) + (8.m14m43 - 3.m24m43 - 8.m13m44 + 3.m23m44)*pow(Z1,2))
572
573	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))
574
575
576	if( prnt )
577	printf "\rDet_b = \r"
578	// 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.
579	printf "TY_qb12 = %15.12g\r",TY_qb12
580	printf "TY_qb21 = %15.12g\r",TY_qb21
581	printf "TY_qb22 = %15.12g\r",TY_qb22
582	printf "TY_qb23 = %15.12g\r",TY_qb23
583	printf "TY_qb32 = %15.12g\r",TY_qb32
584	endif
585
586	// /* Matrix representation of the determinant of the of the system where row refering to
587	// the variable c1 is replaced by solution vector */
588	NVAR TY_qc112 = root:yuk:TY_qc112
589	NVAR TY_qc121 = root:yuk:TY_qc121
590	NVAR TY_qc122 = root:yuk:TY_qc122
591	NVAR TY_qc123 = root:yuk:TY_qc123
592	NVAR TY_qc132 = root:yuk:TY_qc132
593
594	TY_qc112 = -(K1exp(Z1)(9.m24m41 - 9.m14m42 + 3.m11(-12.m14m41 + 4.m24m41 + 8.m14m42 - 3.m24m42) + m44pow(3. - 2.m11,2)))/9.
595
596	TY_qc121 = (K2exp(Z2)(9.m24m31 - 9.m14m32 + 3.m11(-12.m14m31 + 4.m24m31 + 8.m14m32 - 3.m24m32) + m34pow(3. - 2.m11,2)))/9.
597
598	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.
599
600
601	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)))
602	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)))
603	t3 = (8.m42 + 4.m41(-3. + Z2) - 3.m42Z2 + 2.pow(Z2,2))
604	t4 = (6.m44 - 2.m44Z2 + 3.m14*pow(Z2,2))
605	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))
606
607	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.;
608
609
610	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))))
611	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)))
612	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))
613	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)))
614
615	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.;
616
617
618	if( prnt )
619	printf "\rDet_c1 = \r"
620	// 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.
621	printf "TY_qc112 = %15.12g\r",TY_qc112
622	printf "TY_qc121 = %15.12g\r",TY_qc121
623	printf "TY_qc122 = %15.12g\r",TY_qc122
624	printf "TY_qc123 = %15.12g\r",TY_qc123
625	printf "TY_qc132 = %15.12g\r",TY_qc132
626	endif
627
628	// /* Matrix representation of the determinant of the of the system where row refering to
629	// the variable c1 is replaced by solution vector */
630	NVAR TY_qc212 = root:yuk:TY_qc212
631	NVAR TY_qc221 = root:yuk:TY_qc221
632	NVAR TY_qc222 = root:yuk:TY_qc222
633	NVAR TY_qc223 = root:yuk:TY_qc223
634	NVAR TY_qc232 = root:yuk:TY_qc232
635
636	TY_qc212 = (K1exp(Z1)(9m23m41 - 9m13m42 + 3m11(-12m13m41 + 4m23m41 + 8m13m42 - 3m23m42) + m43pow(3 - 2m11,2)))/9.
637
638	TY_qc221 = -(K2exp(Z2)(9m23m31 - 9m13m32 + 3m11(-12m13m31 + 4m23m31 + 8m13m32 - 3m23m32) + m33pow(3 - 2m11,2)))/9.
639
640	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.
641
642
643	t1 = (-(m33(Z1 + Z2)(2m42 + Z2(-2m41 + Z2))) + m32(Z1 + Z2)(2m43 + m13Z2(-2m41 + Z2)) + 2(m13m42 - m43)Z2(-(Z1Z2) + m31*(Z1 + Z2)))
644	t2 = (2(3m41 + 4m11m41 - 3m11m42)Z1pow(Z2,2) + 3m32(Z1 + Z2)(2m41 + m11pow(Z2,2)) - m31(Z1 + Z2)(6m42 + (3 + 4m11)pow(Z2,2)))
645	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))
646
647	TY_qc223 = (2phipow(Z2,-2)(9t1 - 4(-2m43Z1 + m33(Z1 + Z2))pow(m11,2)pow(Z2,2) + 3m23t2 + 6m11t3)*pow(Z1 + Z2,-1))/3.
648
649
650	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))))
651	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)))
652	t3 = (-8m32m43 + m33(m42(8 - 3Z1) + 4m41(-3 + Z1)) - 4m31m43(-3 + Z1) + 3m32m43Z1 - 2(3m13m41 - 2m13m42 + m43)*pow(Z1,2))
653	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)))
654
655	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.
656
657
658	if( prnt )
659	printf "\rDet_c2 = \r"
660	// 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.
661	printf "TY_qc212 = %15.12g\r",TY_qc212
662	printf "TY_qc221 = %15.12g\r",TY_qc221
663	printf "TY_qc222 = %15.12g\r",TY_qc222
664	printf "TY_qc223 = %15.12g\r",TY_qc223
665	printf "TY_qc232 = %15.12g\r",TY_qc232
666	endif
667
668	// /* coefficient matrices of nonlinear equation 1 */
669	NVAR TY_A12 = root:yuk:TY_A12
670	NVAR TY_A21 = root:yuk:TY_A21
671	NVAR TY_A22 = root:yuk:TY_A22
672	NVAR TY_A23 = root:yuk:TY_A23
673	NVAR TY_A32 = root:yuk:TY_A32
674	NVAR TY_A41 = root:yuk:TY_A41
675	NVAR TY_A42 = root:yuk:TY_A42
676	NVAR TY_A43 = root:yuk:TY_A43
677	NVAR TY_A52 = root:yuk:TY_A52
678
679	t1 = (Z1(2TY_qb12(-1 + Z1)(Z1 + Z2) - Z1(2TY_qc212Z1 + TY_qc112(Z1 + Z2))) + TY_qa12(Z1 + Z2)(-2 + pow(Z1,2)))
680	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)))
681
682	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);
683
684
685	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)))
686	t2 = (TY_qb21TY_qc112 + TY_qc121(TY_qa12 + TY_qb12 + 2TY_qc112Z1))
687	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))
688
689	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);
690
691
692	t1 = (TY_qb22TY_qc112 + TY_qc122(TY_qa12 + TY_qb12 + 2TY_qc112Z1))
693	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)))
694	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)))
695
696	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);
697
698
699	t1 = ((TY_q22TY_qc112 + TY_qc123(TY_qc112 + TY_qc212) + TY_qc112TY_qc223)Z1 + TY_qc112TY_qc123Z2)
700	t2 = (TY_qa12TY_qc123 + TY_qa23TY_qc112(1 + Z1) + Z1(TY_qb23TY_qc112 + TY_qc123(TY_qa12 + TY_qb12 + 2TY_qc112Z1)))
701	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)))
702
703	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);
704
705
706	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))
707	t2 = (TY_qc132TY_qc212 + TY_qc123TY_qc221 + TY_qc122TY_qc222 + TY_qc121TY_qc223 + TY_qc112*TY_qc232)
708	t3 = ((TY_q22 + TY_qc132)TY_qc212 + TY_qc123TY_qc221 + TY_qc122TY_qc222 + TY_qc121TY_qc223 + TY_qc112*TY_qc232)
709	t4 = (2TY_qc121TY_qc123 + 2TY_qc112TY_qc132 + pow(TY_qc122,2))
710	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))
711	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))
712	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)
713	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)
714	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))
715	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))
716
717	TY_A32 = exp(-2Z1 - Z2)pow(TY_q22,-2)pow(Z1,-3)t10*pow(Z1 + Z2,-1);
718
719
720	t1 = ((-(TY_qc132TY_qc221) - TY_qc121TY_qc232)Z2 + ((TY_q22 + TY_qc132)TY_qc221 + TY_qc121TY_qc232)(Z1 + Z2)*exp(Z1))
721	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))
722	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))
723	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))
724
725	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);
726
727
728	t1 = (TY_qb32TY_qc122 + (TY_qa22 + TY_qb22)TY_qc132 + TY_qc122(TY_q22 + 2TY_qc132)*Z1)
729	t2 = (TY_qc132TY_qc222 + TY_qc122TY_qc232)
730	t3 = ((TY_q22 + TY_qc132)TY_qc222 + TY_qc122TY_qc232)
731	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)))
732	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)))
733	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)
734
735	TY_A42 = exp(-2Z1 - Z2)pow(TY_q22,-2)pow(Z1,-3)t6*pow(Z1 + Z2,-1);
736
737
738	t1 = ((TY_qc132TY_qc223 + TY_qc123TY_qc232)Z2 - ((TY_q22 + TY_qc132)TY_qc223 + TY_qc123TY_qc232)(Z1 + Z2)*exp(Z1))
739	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))
740	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))
741	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))
742
743	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);
744
745
746	t1 = (TY_qc132Z2 - (TY_q22 + TY_qc132)(Z1 + Z2)*exp(Z1))
747	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)))
748	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)))
749
750	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);
751
752
753	// normalize A
754	// /*double norm_A = sqrt(pow(TY_A52,2)+pow(TY_A43,2)+pow(TY_A42, 2)+pow(TY_A41, 2)+pow(TY_A32, 2)+
755	// pow(TY_A23,2)+pow(TY_A22,2)+pow(TY_A21, 2)+pow(TY_A12, 2));
756	// TY_A12 /= norm_A;
757	// TY_A21 /= norm_A;
758	// TY_A22 /= norm_A;
759	// TY_A23 /= norm_A;
760	// TY_A32 /= norm_A;
761	// TY_A41 /= norm_A;
762	// TY_A42 /= norm_A;
763	// TY_A43 /= norm_A;
764	// TY_A52 /= norm_A;*/
765
766	if( prnt )
767	printf "\rNonlinear equation 1 = \r"
768	// printf "%f\t\t%f\t\t%f\r", 0., TY_A12, 0.
769	// printf "%f\t\t%f\t\t%f\r", TY_A21, TY_A22, TY_A23
770	// printf "%f\t\t%f\t\t%f\r", 0., TY_A32, 0.
771	// printf "%f\t\t%f\t\t%f\r", TY_A41, TY_A42, TY_A43
772	// printf "%f\t\t%f\t\t%f\r", 0., TY_A52, 0.
773	printf "TY_A12 = %15.12g\r",TY_A12
774	printf "TY_A21 = %15.12g\r",TY_A21
775	printf "TY_A22 = %15.12g\r",TY_A22
776	printf "TY_A23 = %15.12g\r",TY_A23
777	printf "TY_A32 = %15.12g\r",TY_A32
778	printf "TY_A41 = %15.12g\r",TY_A41
779	printf "TY_A42 = %15.12g\r",TY_A42
780	printf "TY_A43 = %15.12g\r",TY_A43
781	printf "TY_A52 = %15.12g\r",TY_A52
782	endif
783
784	// /* coefficient matrices of nonlinear equation 2 */
785	NVAR TY_B12 = root:yuk:TY_B12
786	NVAR TY_B14 = root:yuk:TY_B14
787	NVAR TY_B21 = root:yuk:TY_B21
788	NVAR TY_B22 = root:yuk:TY_B22
789	NVAR TY_B23 = root:yuk:TY_B23
790	NVAR TY_B24 = root:yuk:TY_B24
791	NVAR TY_B25 = root:yuk:TY_B25
792	NVAR TY_B32 = root:yuk:TY_B32
793	NVAR TY_B34 = root:yuk:TY_B34
794
795
796
797	t1 = (TY_qa12TY_qc221 + TY_qa21TY_qc212(1 + Z2) + Z2(TY_qb21TY_qc212 + TY_qc221(TY_qa12 + TY_qb12 + 2TY_qc212Z2)))
798	t2 = (-(TY_qc121TY_qc212) - TY_qc112TY_qc221)
799	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))
800	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))
801
802	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);
803
804
805
806
807	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))
808	t2 = (TY_qc123TY_qc212 + TY_qc112(TY_q22 + TY_qc223))
809	t3 = (TY_qa23TY_qc212 + TY_qa12(TY_q22 + TY_qc223))
810	t4 = (TY_qa23TY_qc212 + TY_qa12(TY_q22 + TY_qc223) + (TY_qb23TY_qc212 + TY_qb12(TY_q22 + TY_qc223))*Z1)
811	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)
812	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)
813
814	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);
815
816
817
818	t1 = (TY_qc221(Z1 + Z2)exp(Z1)*pow(Z2,2))
819	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))
820	t3 = (-(TY_qc121Z1pow(Z2,2)) + (Z1 + Z2)exp(Z1)(TY_qa21 + (TY_qa21 + TY_qb21)Z2 + TY_qc221pow(Z2,2)))
821
822	TY_B21 = 6phiTY_qc221exp(-Z1 - 2Z2)pow(TY_q22,-2)pow(Z2,-3)(-t1 + exp(2Z2)t2 + 2exp(Z2)t3)pow(Z1 + Z2,-1);
823
824
825
826	t1 = (TY_qb22TY_qc221 + TY_qc222(TY_qa21 + TY_qb21 + 2TY_qc221Z2))
827	t2 = ((Z1 + Z2)(TY_qa21TY_qc222 + TY_qa22TY_qc221(1 + Z2) + Z2t1)exp(Z1) + (-(TY_qc122TY_qc221) - TY_qc121TY_qc222)Z1pow(Z2,2))
828	t3 = (TY_qc122TY_qc221 + TY_qc121TY_qc222)
829	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))
830	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))
831
832	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);
833
834
835
836
837	t1 = (2TY_qc221TY_qc223 + 2TY_qc212TY_qc232 + pow(TY_qc222,2))
838	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)))
839	t3 = (-(TY_qc132TY_qc212) - TY_qc123TY_qc221 - TY_qc122TY_qc222 - TY_qc121TY_qc223 - TY_qc112*TY_qc232)
840	t4 = (TY_qc132TY_qc212 + TY_qc123TY_qc221 + TY_qc122TY_qc222 + TY_qc121(TY_q22 + TY_qc223) + TY_qc112*TY_qc232)
841	t5 = (TY_qa32TY_qc212 + TY_qa23TY_qc221 + TY_qa22TY_qc222 + TY_qa21(TY_q22 + TY_qc223) + TY_qa12*TY_qc232)
842	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)
843	t7 = (TY_qb32TY_qc212 + TY_qb23TY_qc221 + TY_qb22TY_qc222 + TY_qb21(TY_q22 + TY_qc223) + TY_qb12*TY_qc232)
844	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)
845	t9 = (TY_qa21 + 2TY_qb21 - 2TY_qc121 - 2TY_qc212 - 2TY_qc221)
846	t10 = (TY_qa22 + 2TY_qb22 - 2TY_qc122 - TY_qc222)
847	t11 = (TY_qa23 + 2TY_qb23 - 2TY_qc123 - 2*TY_qc223)
848	t12 = (TY_qa32 + 2TY_qb32 - 2TY_qc132 - 2*TY_qc232)
849	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))
850	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)
851	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)))))
852
853	TY_B23 = exp(-Z1 - 2Z2)pow(TY_q22,-2)pow(Z2,-3)t15*pow(Z1 + Z2,-1);
854
855
856	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))
857	t2 = (TY_qc123TY_qc222 + TY_qc122TY_qc223)
858	t3 = (TY_qc123TY_qc222 + TY_qc122(TY_q22 + TY_qc223))
859	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)))
860	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)))
861	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)
862
863	TY_B24 = exp(-Z1 - 2Z2)pow(TY_q22,-2)pow(Z2,-3)t6*pow(Z1 + Z2,-1);
864
865
866	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))
867	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))))
868
869	TY_B25 = -6phiexp(-Z1 - 2Z2)pow(TY_q22,-2)pow(Z2,-3)t2*pow(Z1 + Z2,-1);
870
871
872	t1 = (TY_qa21TY_qc232 + TY_qa32TY_qc221(1 + Z2) + Z2(TY_qb32TY_qc221 + TY_qc232(TY_qa21 + TY_qb21 + 2TY_qc221Z2)))
873	t2 = (-(TY_qc132TY_qc221) - TY_qc121TY_qc232)
874	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))
875	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))
876
877	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);
878
879
880	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))
881	t2 = (TY_qc132(TY_q22 + TY_qc223) + TY_qc123TY_qc232)
882	t3 = (TY_qa32(TY_q22 + TY_qc223) + TY_qa23TY_qc232)
883	t4 = (TY_qa32(TY_q22 + TY_qc223) + TY_qa23TY_qc232 + (TY_qb32(TY_q22 + TY_qc223) + TY_qb23TY_qc232)*Z1)
884	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)
885	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))
886
887	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);
888
889
890	// /*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));
891	//
892	// TY_B12 /= norm_B;
893	// TY_B14 /= norm_B;
894	// TY_B21 /= norm_B;
895	// TY_B22 /= norm_B;
896	// TY_B23 /= norm_B;
897	// TY_B24 /= norm_B;
898	// TY_B25 /= norm_B;
899	// TY_B32 /= norm_B;
900	// TY_B34 /= norm_B; */
901
902	if( prnt )
903	printf "\rNonlinear equation 2 = \r"
904	// printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0., TY_B12, 0., TY_B14, 0.
905	// 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
906	// printf "%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0., TY_B32, 0., TY_B34, 0.
907	printf "TY_B12 = %15.12g\r",TY_B12
908	printf "TY_B14 = %15.12g\r",TY_B14
909	printf "TY_B21 = %15.12g\r",TY_B21
910	printf "TY_B22 = %15.12g\r",TY_B22
911	printf "TY_B23 = %15.12g\r",TY_B23
912	printf "TY_B24 = %15.12g\r",TY_B24
913	printf "TY_B25 = %15.12g\r",TY_B25
914	printf "TY_B32 = %15.12g\r",TY_B32
915	printf "TY_B34 = %15.12g\r",TY_B34
916	endif
917
918	// /* decrease order of nonlinear equation 1 by means of equation 2 */
919	NVAR TY_F14 = root:yuk:TY_F14
920	NVAR TY_F16 = root:yuk:TY_F16
921	NVAR TY_F18 = root:yuk:TY_F18
922	NVAR TY_F23 = root:yuk:TY_F23
923	NVAR TY_F24 = root:yuk:TY_F24
924	NVAR TY_F25 = root:yuk:TY_F25
925	NVAR TY_F26 = root:yuk:TY_F26
926	NVAR TY_F27 = root:yuk:TY_F27
927	NVAR TY_F28 = root:yuk:TY_F28
928	NVAR TY_F29 = root:yuk:TY_F29
929	NVAR TY_F32 = root:yuk:TY_F32
930	NVAR TY_F33 = root:yuk:TY_F33
931	NVAR TY_F34 = root:yuk:TY_F34
932	NVAR TY_F35 = root:yuk:TY_F35
933	NVAR TY_F36 = root:yuk:TY_F36
934	NVAR TY_F37 = root:yuk:TY_F37
935	NVAR TY_F38 = root:yuk:TY_F38
936	NVAR TY_F39 = root:yuk:TY_F39
937	NVAR TY_F310 = root:yuk:TY_F310
938
939	TY_F14 = -(TY_A32TY_B12TY_B32) + TY_A52pow(TY_B12,2) + TY_A12pow(TY_B32,2)
940	TY_F16 = 2TY_A52TY_B12TY_B14 - TY_A32TY_B14TY_B32 - TY_A32TY_B12TY_B34 + 2TY_A12TY_B32TY_B34
941	TY_F18 = -(TY_A32TY_B14TY_B34) + TY_A52pow(TY_B14,2) + TY_A12pow(TY_B34,2)
942	TY_F23 = 2TY_A52TY_B12TY_B21 - TY_A41TY_B12TY_B32 - TY_A32TY_B21TY_B32 + TY_A21pow(TY_B32,2)
943	TY_F24 = 2TY_A52TY_B12TY_B22 - TY_A42TY_B12TY_B32 - TY_A32TY_B22TY_B32 + TY_A22pow(TY_B32,2)
944	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)
945	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
946	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)
947	TY_F28 = 2TY_A52TY_B14TY_B24 - TY_A42TY_B14TY_B34 - TY_A32TY_B24TY_B34 + TY_A22pow(TY_B34,2)
948	TY_F29 = 2TY_A52TY_B14TY_B25 - TY_A43TY_B14TY_B34 - TY_A32TY_B25TY_B34 + TY_A23pow(TY_B34,2)
949	TY_F32 = -(TY_A41TY_B21TY_B32) + TY_A52*pow(TY_B21,2)
950	TY_F33 = 2TY_A52TY_B21TY_B22 - TY_A42TY_B21TY_B32 - TY_A41TY_B22*TY_B32
951	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)
952	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
953	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)
954	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
955	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)
956	TY_F39 = 2TY_A52TY_B24TY_B25 - TY_A43TY_B24TY_B34 - TY_A42TY_B25*TY_B34
957	TY_F310 = -(TY_A43TY_B25TY_B34) + TY_A52*pow(TY_B25,2)
958
959	if( prnt )
960	printf "\rF = \r"
961	// 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.
962	// 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.
963	// 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
964	printf "TY_F14 = %15.12g\r",TY_F14
965	printf "TY_F16 = %15.12g\r",TY_F16
966	printf "TY_F18 = %15.12g\r",TY_F18
967	printf "TY_F23 = %15.12g\r",TY_F23
968	printf "TY_F24 = %15.12g\r",TY_F24
969	printf "TY_F25 = %15.12g\r",TY_F25
970	printf "TY_F26 = %15.12g\r",TY_F26
971	printf "TY_F27 = %15.12g\r",TY_F27
972	printf "TY_F28 = %15.12g\r",TY_F28
973	printf "TY_F29 = %15.12g\r",TY_F29
974	printf "TY_F32 = %15.12g\r",TY_F32
975	printf "TY_F33 = %15.12g\r",TY_F33
976	printf "TY_F34 = %15.12g\r",TY_F34
977	printf "TY_F35 = %15.12g\r",TY_F35
978	printf "TY_F36 = %15.12g\r",TY_F36
979	printf "TY_F37 = %15.12g\r",TY_F37
980	printf "TY_F38 = %15.12g\r",TY_F38
981	printf "TY_F39 = %15.12g\r",TY_F39
982	printf "TY_F310 = %15.12g\r",TY_F310
983	endif
984
985	NVAR TY_G13 = root:yuk:TY_G13
986	NVAR TY_G14 = root:yuk:TY_G14
987	NVAR TY_G15 = root:yuk:TY_G15
988	NVAR TY_G16 = root:yuk:TY_G16
989	NVAR TY_G17 = root:yuk:TY_G17
990	NVAR TY_G18 = root:yuk:TY_G18
991	NVAR TY_G19 = root:yuk:TY_G19
992	NVAR TY_G110 = root:yuk:TY_G110
993	NVAR TY_G111 = root:yuk:TY_G111
994	NVAR TY_G112 = root:yuk:TY_G112
995	NVAR TY_G113 = root:yuk:TY_G113
996	NVAR TY_G22 = root:yuk:TY_G22
997	NVAR TY_G23 = root:yuk:TY_G23
998	NVAR TY_G24 = root:yuk:TY_G24
999	NVAR TY_G25 = root:yuk:TY_G25
1000	NVAR TY_G26 = root:yuk:TY_G26
1001	NVAR TY_G27 = root:yuk:TY_G27
1002	NVAR TY_G28 = root:yuk:TY_G28
1003	NVAR TY_G29 = root:yuk:TY_G29
1004	NVAR TY_G210 = root:yuk:TY_G210
1005	NVAR TY_G211 = root:yuk:TY_G211
1006	NVAR TY_G212 = root:yuk:TY_G212
1007	NVAR TY_G213 = root:yuk:TY_G213
1008	NVAR TY_G214 = root:yuk:TY_G214
1009
1010
1011	TY_G13 = -(TY_B12*TY_F32)
1012	TY_G14 = -(TY_B12*TY_F33)
1013	TY_G15 = TY_B32TY_F14 - TY_B14TY_F32 - TY_B12*TY_F34
1014	TY_G16 = -(TY_B14TY_F33) - TY_B12TY_F35
1015	TY_G17 = TY_B34TY_F14 + TY_B32TY_F16 - TY_B14TY_F34 - TY_B12TY_F36
1016	TY_G18 = -(TY_B14TY_F35) - TY_B12TY_F37
1017	TY_G19 = TY_B34TY_F16 + TY_B32TY_F18 - TY_B14TY_F36 - TY_B12TY_F38
1018	TY_G110 = -(TY_B14TY_F37) - TY_B12TY_F39
1019	TY_G111 = TY_B34TY_F18 - TY_B12TY_F310 - TY_B14*TY_F38
1020	TY_G112 = -(TY_B14*TY_F39)
1021	TY_G113 = -(TY_B14*TY_F310)
1022	TY_G22 = -(TY_B21*TY_F32)
1023	TY_G23 = -(TY_B22TY_F32) - TY_B21TY_F33
1024	TY_G24 = TY_B32TY_F23 - TY_B23TY_F32 - TY_B22TY_F33 - TY_B21TY_F34
1025	TY_G25 = TY_B32TY_F24 - TY_B24TY_F32 - TY_B23TY_F33 - TY_B22TY_F34 - TY_B21*TY_F35
1026	TY_G26 = TY_B34TY_F23 + TY_B32TY_F25 - TY_B25TY_F32 - TY_B24TY_F33 - TY_B23TY_F34 - TY_B22TY_F35 - TY_B21*TY_F36
1027	TY_G27 = TY_B34TY_F24 + TY_B32TY_F26 - TY_B25TY_F33 - TY_B24TY_F34 - TY_B23TY_F35 - TY_B22TY_F36 - TY_B21*TY_F37
1028	TY_G28 = TY_B34TY_F25 + TY_B32TY_F27 - TY_B25TY_F34 - TY_B24TY_F35 - TY_B23TY_F36 - TY_B22TY_F37 - TY_B21*TY_F38
1029	TY_G29 = TY_B34TY_F26 + TY_B32TY_F28 - TY_B25TY_F35 - TY_B24TY_F36 - TY_B23TY_F37 - TY_B22TY_F38 - TY_B21*TY_F39
1030	TY_G210 = TY_B34TY_F27 + TY_B32TY_F29 - TY_B21TY_F310 - TY_B25TY_F36 - TY_B24TY_F37 - TY_B23TY_F38 - TY_B22*TY_F39
1031	TY_G211 = TY_B34TY_F28 - TY_B22TY_F310 - TY_B25TY_F37 - TY_B24TY_F38 - TY_B23*TY_F39
1032	TY_G212 = TY_B34TY_F29 - TY_B23TY_F310 - TY_B25TY_F38 - TY_B24TY_F39
1033	TY_G213 = -(TY_B24TY_F310) - TY_B25TY_F39
1034	TY_G214 = -(TY_B25*TY_F310)
1035
1036	if( prnt )
1037	printf "\rG = \r"
1038	// 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.
1039	// 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
1040	printf "TY_G13 = %15.12g\r",TY_G13
1041	printf "TY_G14 = %15.12g\r",TY_G14
1042	printf "TY_G15 = %15.12g\r",TY_G15
1043	printf "TY_G16 = %15.12g\r",TY_G16
1044	printf "TY_G17 = %15.12g\r",TY_G17
1045	printf "TY_G18 = %15.12g\r",TY_G18
1046	printf "TY_G19 = %15.12g\r",TY_G19
1047	printf "TY_G110 = %15.12g\r",TY_G110
1048	printf "TY_G111 = %15.12g\r",TY_G111
1049	printf "TY_G112 = %15.12g\r",TY_G112
1050	printf "TY_G113 = %15.12g\r",TY_G113
1051	printf "TY_G22 = %15.12g\r",TY_G22
1052	printf "TY_G23 = %15.12g\r",TY_G23
1053	printf "TY_G24 = %15.12g\r",TY_G24
1054	printf "TY_G25 = %15.12g\r",TY_G25
1055	printf "TY_G26 = %15.12g\r",TY_G26
1056	printf "TY_G27 = %15.12g\r",TY_G27
1057	printf "TY_G28 = %15.12g\r",TY_G28
1058	printf "TY_G29 = %15.12g\r",TY_G29
1059	printf "TY_G210 = %15.12g\r",TY_G210
1060	printf "TY_G211 = %15.12g\r",TY_G211
1061	printf "TY_G212 = %15.12g\r",TY_G212
1062	printf "TY_G213 = %15.12g\r",TY_G213
1063	printf "TY_G214 = %15.12g\r",TY_G214
1064	endif
1065
1066	Make/O/D/N=23 TY_w
1067
1068	// coefficients for polynomial
1069	TY_w[0] = (-(TY_A21TY_B12) + TY_A12TY_B21)(TY_A52TY_B21 - TY_A41TY_B32)pow(TY_B21,2)*pow(TY_B32,3)
1070
1071	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
1072	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
1073
1074	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
1075	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
1076	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)
1077
1078	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
1079	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
1080	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
1081	TY_w[3] += 2TY_B12TY_G22*TY_G27;
1082
1083	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
1084	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
1085	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
1086	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)
1087	TY_w[4] += TY_B12*pow(TY_G25,2)
1088
1089	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
1090	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
1091	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
1092	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
1093	TY_w[5] += -TY_B22TY_G13TY_G28 - TY_B21TY_G14TY_G28 + 2TY_B12TY_G23TY_G28 - TY_B21TY_G13TY_G29 + 2TY_B12TY_G22TY_G29
1094
1095	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
1096	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
1097	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
1098	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
1099	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))
1100	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)
1101
1102	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
1103	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
1104	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
1105	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
1106	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
1107	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
1108	TY_w[7] += 2TY_B12TY_G24*TY_G29
1109
1110	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
1111	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
1112	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
1113	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
1114	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
1115	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
1116	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)
1117	TY_w[8] += TY_B12*pow(TY_G27,2)
1118
1119	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
1120	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
1121	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
1122	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
1123	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
1124	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
1125	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
1126	TY_w[9] += -TY_B23TY_G15TY_G29 - TY_B22TY_G16TY_G29 - TY_B21TY_G17TY_G29 + 2TY_B14TY_G24TY_G29 + 2TY_B12TY_G26TY_G29
1127
1128	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
1129	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
1130	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
1131	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
1132	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
1133	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
1134	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
1135	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))
1136	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)
1137
1138	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
1139	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
1140	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
1141	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
1142	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
1143	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
1144	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
1145	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
1146	TY_w[11] += 2TY_B14TY_G26TY_G29 + 2TY_B12TY_G28TY_G29
1147
1148	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
1149	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
1150	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
1151	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
1152	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
1153	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
1154	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
1155	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))
1156	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)
1157
1158	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
1159	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
1160	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
1161	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
1162	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
1163	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
1164	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
1165	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
1166
1167	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
1168	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
1169	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
1170	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
1171	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
1172	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
1173	TY_w[14] += 2TY_B12TY_G211TY_G29 + TY_B32(2(TY_G113TY_G17 + TY_G112TY_G18 + TY_G111TY_G19) + pow(TY_G110,2))
1174	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)
1175
1176	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
1177	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
1178	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
1179	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
1180	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
1181	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
1182	TY_w[15] += -TY_B21TY_G113TY_G29 - TY_B25TY_G19TY_G29 + 2TY_B14TY_G210TY_G29 + 2TY_B12TY_G212TY_G29
1183
1184	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
1185	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
1186	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
1187	TY_w[16] += 2TY_B14TY_G213TY_G27 - TY_B25TY_G111TY_G28 - TY_B24TY_G112TY_G28 - TY_B23TY_G113TY_G28 + 2TY_B14TY_G212TY_G28
1188	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
1189	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))
1190	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)
1191
1192	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
1193	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
1194	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
1195	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
1196	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)
1197
1198	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
1199	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
1200	TY_w[18] += -TY_B25TY_G18TY_G213 - TY_B24TY_G19TY_G213 + 2TY_B12TY_G211TY_G213 - TY_B25TY_G113*TY_G28
1201	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
1202	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)
1203
1204	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
1205	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
1206	TY_w[19] += -TY_B21TY_G113TY_G213 - TY_B25TY_G19TY_G213 + 2TY_B14TY_G210TY_G213 + 2TY_B12TY_G212TY_G213 - TY_B25TY_G113TY_G29
1207	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)
1208
1209	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
1210	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
1211	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)
1212
1213	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)
1214
1215	TY_w[22] = (-(TY_A23TY_B14) + TY_A12TY_B25)(TY_A52TY_B25 - TY_A43TY_B34)pow(TY_B25,2)*pow(TY_B34,3)
1216
1217	if( prnt )
1218	printf "\rCoefficients of polynomial\r"
1219	variable i
1220	for ( i = 0; i < 23; i+=1 )
1221	printf "w[%d] = %g\r", i, TY_w[i]
1222	endfor
1223	printf "\r"
1224	endif
1225
1226	end
1227
1228	Function TY_capQ( d2 )
1229	Variable d2
1230
1231	NVAR TY_B32 = root:yuk:TY_B32
1232	NVAR TY_B34 = root:yuk:TY_B34
1233
1234	return d2 * TY_B32 + pow( d2, 3 ) * TY_B34
1235	end
1236
1237	Function TY_V( d2 )
1238	variable d2
1239
1240	NVAR TY_G13 = root:yuk:TY_G13
1241	NVAR TY_G14 = root:yuk:TY_G14
1242	NVAR TY_G15 = root:yuk:TY_G15
1243	NVAR TY_G16 = root:yuk:TY_G16
1244	NVAR TY_G17 = root:yuk:TY_G17
1245	NVAR TY_G18 = root:yuk:TY_G18
1246	NVAR TY_G19 = root:yuk:TY_G19
1247	NVAR TY_G110 = root:yuk:TY_G110
1248	NVAR TY_G111 = root:yuk:TY_G111
1249	NVAR TY_G112 = root:yuk:TY_G112
1250	NVAR TY_G113 = root:yuk:TY_G113
1251
1252	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 )
1253	end
1254
1255	Function TY_capW( d2 )
1256	Variable d2
1257
1258	variable tmp
1259
1260	NVAR TY_G22 = root:yuk:TY_G22
1261	NVAR TY_G23 = root:yuk:TY_G23
1262	NVAR TY_G24 = root:yuk:TY_G24
1263	NVAR TY_G25 = root:yuk:TY_G25
1264	NVAR TY_G26 = root:yuk:TY_G26
1265	NVAR TY_G27 = root:yuk:TY_G27
1266	NVAR TY_G28 = root:yuk:TY_G28
1267	NVAR TY_G29 = root:yuk:TY_G29
1268	NVAR TY_G210 = root:yuk:TY_G210
1269	NVAR TY_G211 = root:yuk:TY_G211
1270	NVAR TY_G212 = root:yuk:TY_G212
1271	NVAR TY_G213 = root:yuk:TY_G213
1272	NVAR TY_G214 = root:yuk:TY_G214
1273
1274
1275	tmp = d2 * TY_G22 + pow( d2, 2 ) * TY_G23 + pow( d2, 3 ) * TY_G24 + pow( d2, 4 ) * TY_G25 + pow( d2, 5 ) * TY_G26
1276	tmp += pow( d2, 6 ) * TY_G27 + pow( d2, 7 ) * TY_G28 + pow( d2, 8 ) * TY_G29 + pow( d2, 9 ) * TY_G210
1277	tmp += pow( d2, 10 ) * TY_G211 + pow( d2, 11 ) * TY_G212 + pow( d2, 12 ) * TY_G213 + pow( d2, 13 ) * TY_G214
1278
1279	return tmp
1280	end
1281
1282	Function TY_X( d2 )
1283	Variable d2
1284
1285	return TY_V( d2 ) / TY_capW( d2 )
1286	end
1287
1288	// solve the linear system depending on d1, d2 using Cramer's rule
1289	//
1290	// a,b,c1,c2 are passed by reference and returned
1291	//
1292	Function TY_SolveLinearEquations( d1, d2, a, b, c1, c2)
1293	Variable d1, d2, &a, &b, &c1, &c2
1294
1295
1296	NVAR TY_q22 = root:yuk:TY_q22
1297	NVAR TY_qa12 = root:yuk:TY_qa12
1298	NVAR TY_qa21 = root:yuk:TY_qa21
1299	NVAR TY_qa22 = root:yuk:TY_qa22
1300	NVAR TY_qa23 = root:yuk:TY_qa23
1301	NVAR TY_qa32 = root:yuk:TY_qa32
1302
1303	NVAR TY_qb12 = root:yuk:TY_qb12
1304	NVAR TY_qb21 = root:yuk:TY_qb21
1305	NVAR TY_qb22 = root:yuk:TY_qb22
1306	NVAR TY_qb23 = root:yuk:TY_qb23
1307	NVAR TY_qb32 = root:yuk:TY_qb32
1308
1309	NVAR TY_qc112 = root:yuk:TY_qc112
1310	NVAR TY_qc121 = root:yuk:TY_qc121
1311	NVAR TY_qc122 = root:yuk:TY_qc122
1312	NVAR TY_qc123 = root:yuk:TY_qc123
1313	NVAR TY_qc132 = root:yuk:TY_qc132
1314
1315	NVAR TY_qc212 = root:yuk:TY_qc212
1316	NVAR TY_qc221 = root:yuk:TY_qc221
1317	NVAR TY_qc222 = root:yuk:TY_qc222
1318	NVAR TY_qc223 = root:yuk:TY_qc223
1319	NVAR TY_qc232 = root:yuk:TY_qc232
1320
1321	Variable det = TY_q22 * d1 * d2
1322	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
1323	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
1324	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
1325	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
1326
1327	a = det_a / det
1328	b = det_b / det
1329	c1 = det_c1 / det
1330	c2 = det_c2 / det
1331	end
1332
1333	//Solve the system of linear and nonlinear equations for given Zi, Ki, phi which gives at
1334	// most 22 solutions for the parameters a,b,ci,di. From the set of solutions choose the
1335	// physical one and return it.
1336	//
1337	//
1338	// a,b,c1,c2,d1,d2 are passed by reference and returned
1339	//
1340	Function TY_SolveEquations( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2, prnt )
1341	Variable Z1, Z2, K1, K2, phi, &a, &b, &c1, &c2, &d1, &d2, prnt
1342
1343
1344	// reduce system to a polynomial from which all solution are extracted
1345	// by doing that a lot of global background variables are set
1346	TY_ReduceNonlinearSystem( Z1, Z2, K1, K2, phi, prnt )
1347
1348	// the two coupled non-linear eqautions were reduced to a
1349	// 22nd order polynomial, the roots are give all possible solutions
1350	// for d2, than d1 can be computed by the function X
1351
1352	Make/O/D/N=23 real_coefficient,imag_coefficient
1353	Make/O/D/N=22 real_root,imag_root
1354
1355	//integer degree of polynomial
1356	variable degree = 22
1357	Variable i
1358
1359	WAVE TY_w = TY_w
1360
1361
1362	////
1363	// now I need to replace this solution with FindRoots/P to get the polynomial roots
1364	////
1365
1366	// vector of real and imaginary coefficients in order of INCREASING powers
1367	for ( i = 0; i <= degree; i+=1 )
1368	// the global variablw TY_w was set by TY_ReduceNonlinearSystem
1369	real_coefficient[i] = TY_w[i]
1370	// imag_coefficient[i] = 0.;
1371	endfor
1372
1373	// zrhqr(real_coefficient, degree, NR_r, NR_i);
1374
1375	FindRoots/P=real_coefficient
1376
1377	WAVE/C W_polyRoots = W_polyRoots
1378
1379	for(i=0; i<degree; i+=1)
1380	real_root[i] = real(W_polyRoots[i])
1381	imag_root[i] = imag(W_polyRoots[i])
1382	endfor
1383
1384	//end - NR solution of polynomial
1385
1386
1387	// show the result if in debug mode
1388	Variable x, y
1389	if ( prnt )
1390	for ( i = 0; i < degree; i+=1 )
1391	x = real_root[i]
1392	y = imag_root[i]
1393	if ( chop( y ) == 0 )
1394	printf "root(%d) = %g\r", i+1, x
1395	else
1396	printf "root(%d) = %g + %g i\r", i+1, x, y
1397	endif
1398	endfor
1399	printf "\r"
1400	endif
1401
1402
1403
1404	// select real roots and those satisfying Q(x) != 0 and W(x) != 0
1405	// Paper: Cluster formation in two-Yukawa Fluids, J. Chem. Phys. 122, 2005
1406	// The right set of (a, b, c1, c2, d1, d2) should have the following properties:
1407	// (1) a > 0
1408	// (2) d1, d2 are real
1409	// (3) vi/Ki > 0 <=> g(Zi) > 0
1410	// (4) if there is still more than root, calculate g(r) for each root
1411	// and g(r) of the correct root should have the minimum average value
1412	// inside the hardcore
1413	Variable var_a, var_b, var_c1, var_c2, var_d1, var_d2
1414	Make/O/D/N=22 sol_a, sol_b, sol_c1, sol_c2, sol_d1, sol_d2
1415
1416	Variable j = 0
1417	for ( i = 0; i < degree; i+=1 )
1418
1419	x = real_root[i]
1420	y = imag_root[i]
1421
1422	if ( chop( y ) == 0 && TY_capW( x ) != 0 && TY_capQ( x ) != 0 )
1423
1424	var_d1 = TY_X( x )
1425	var_d2 = x
1426
1427	// solution of linear system for given d1, d2 to obtain a,b,ci,di
1428	// var_a, var_b, var_c1, var_c2 passed by reference
1429	TY_SolveLinearEquations( var_d1, var_d2, var_a, var_b, var_c1, var_c2 )
1430
1431	// select physical solutions, for details check paper: "Cluster formation in
1432	// two-Yukawa fluids", J. Chem. Phys. 122 (2005)
1433	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 )
1434	sol_a[j] = var_a
1435	sol_b[j] = var_b
1436	sol_c1[j] = var_c1
1437	sol_c2[j] = var_c2
1438	sol_d1[j] = var_d1
1439	sol_d2[j] = var_d2
1440
1441	if ( prnt )
1442	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] ) )
1443	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] ) )
1444	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] ) )
1445	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] ) )
1446	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] ) )
1447	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] ) )
1448
1449	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
1450	endif
1451
1452	j+=1
1453	endif //var_a >0...
1454	endif //chop
1455	endfor
1456	// number remaining roots
1457	Variable n_roots = j
1458
1459	// if there is still more than one root left, than choose the one with the minimum
1460	// average value inside the hardcore
1461	if ( n_roots > 1 )
1462
1463	/////
1464	// it seems like this section should all be replaced in bulk with internal FFT code, rather than slow integration
1465	//
1466	// -- also, be sure to handle r=0, or the sum will always be INF
1467	////
1468
1469	// the number of q values should be a power of 2
1470	// in order to speed up the FFT
1471	/// int n = 1 << 14;
1472	Variable n=16384 //2^14 points
1473
1474	// the maximum q value should be large enough
1475	// to enable a reasoble approximation of g(r)
1476	variable qmax = 16 * 10 * 2 * pi
1477	Variable q, dq = qmax / ( n - 1 )
1478
1479	// step size for g(r)
1480	variable dr
1481
1482	// allocate memory for pair correlation function g(r)
1483	// and structure factor S(q)
1484	Make/O/D/N=(n) sq,gr //gr will be redimensioned!!
1485
1486	// loop over all remaining roots
1487	Variable minVal = 1e50 //a really big number
1488	Variable selected_root = 10
1489	Variable sumVal = 0
1490
1491	for ( j = 0; j < n_roots; j+=1)
1492
1493	// calculate structure factor at different q values
1494	for ( i = 0; i < n; i+=1)
1495
1496	q = dq * i
1497	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] )
1498
1499	if(i<10 && prnt)
1500	printf "after SqTwoYukawa: s(q) = %g\r",sq[i]
1501	endif
1502
1503	endfor
1504
1505	// calculate pair correlation function for given
1506	// structure factor, g(r) is computed at values
1507	// r(i) = i * dr
1508
1509	// Yuk_SqToGr( phi, dq, sq, dr, gr, n )
1510
1511
1512	Yuk_SqToGr_FFT( phi, dq, sq, dr, gr, n )
1513
1514	// determine sum inside the hardcore
1515	// 0 =< r < 1 of the pair-correlation function
1516	sumVal = 0
1517	for (i = 0; i < floor( 1. / dr ); i+=1 )
1518
1519	sumVal += abs( gr[i] )
1520
1521	if(i<10 && prnt)
1522	printf "g(r) in core = %g\r",abs(gr[i])
1523	endif
1524
1525	endfor
1526
1527	if ( sumVal < minVal )
1528	minVal = sumVal
1529	selected_root = j
1530	endif
1531
1532	if(prnt)
1533	printf "min = %g sum = %g\r",minVal,sumVal
1534	endif
1535
1536	endfor
1537
1538
1539	// physical solution was found
1540	a = sol_a [ selected_root ] //sol_a [ selected_root ];
1541	b = sol_b [ selected_root ]
1542	c1 = sol_c1[ selected_root ]
1543	c2 = sol_c2[ selected_root ]
1544	d1 = sol_d1[ selected_root ]
1545	d2 = sol_d2[ selected_root ]
1546
1547	return 1
1548
1549	else
1550	if ( n_roots == 1 )
1551
1552	a = sol_a [0]
1553	b = sol_b [0]
1554	c1 = sol_c1[0]
1555	c2 = sol_c2[0]
1556	d1 = sol_d1[0]
1557	d2 = sol_d2[0]
1558
1559	return 1
1560	else
1561	// no solution was found
1562	return 0
1563	endif
1564	endif
1565
1566	end
1567
1568
1569
1570
1571	//
1572	Function Yuk_SqToGr_FFT( phi, dq, sq, dr, gr, n )
1573	Variable phi, dq
1574	WAVE sq
1575	Variable &dr
1576	WAVE gr
1577	Variable n
1578
1579	Variable npts,ii,rval,jj,qval,spread=1
1580	Variable alpha
1581
1582
1583	WaveStats/Q sq
1584	npts = V_npnts
1585
1586	dr = 2pi/(nptsdq)
1587
1588	Make/O/D/N=(npts) temp
1589
1590	temp = p*(sq[p] - 1)
1591	alpha = npts * pow( dq, 3 ) / ( 24 * pi * pi * phi )
1592
1593	FFT/OUT=1/DEST=W_FFT temp
1594
1595
1596	WAVE/C W_FFT = W_FFT
1597
1598	Redimension/N=(numpnts(W_FFT)) gr
1599
1600	gr = 1 + alpha/p*imag(W_FFT)
1601
1602	gr[0] = 0
1603
1604	// Killwaves/Z temp
1605
1606	return(0)
1607
1608	End
1609
1610	////////////////////////end converted procedures //////////////////////////////////

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

Download in other formats:

Original Format