source: sasview/sansmodels/src/libigor/2Y_TwoYukawa.c @ 8037f4b

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 8037f4b was a42ea15, checked in by Jae Cho <jhjcho@…>, 12 years ago

Yukawa structure model kernels (IGOR/NIST)
  • Property mode set to 100644
File size: 73.3 KB
Line 
1/*
2 *  twoyukawa.c
3 *  twoyukawa
4 *
5 *  Created by Marcus Hennig on 5/7/10.
6 *  Copyright 2010 __MyCompanyName__. All rights reserved.
7 *
8 */
9
10#include "2Y_TwoYukawa.h"
11#include "2Y_cpoly.h"
12#include "2Y_utility.h"
13#include "2Y_PairCorrelation.h"
14#include <stdio.h>
15#include <math.h>
16#include <stdlib.h>
17
18/*
19 ==================================================================================================
20 
21 The two-yukawa structure factor is uniquley determined by 6 parameters a, b, c1, c2, d1, d2,
22 which are the solution of a system of 6 equations ( 4 linear, 2 nonlinear ). The solution can
23 constructed by the roots of a polynomial of 22nd degree. For more details see attached
24 Mathematica snotebook, where a derivation is given
25 
26 ==================================================================================================
27 */
28
29double TY_q22;
30double TY_qa12, TY_qa21, TY_qa22, TY_qa23, TY_qa32;
31double TY_qb12, TY_qb21, TY_qb22, TY_qb23, TY_qb32;
32double TY_qc112, TY_qc121, TY_qc122, TY_qc123, TY_qc132;
33double TY_qc212, TY_qc221, TY_qc222, TY_qc223, TY_qc232;
34double TY_A12, TY_A21, TY_A22, TY_A23, TY_A32, TY_A41, TY_A42, TY_A43, TY_A52;
35double TY_B12, TY_B14, TY_B21, TY_B22, TY_B23, TY_B24, TY_B25, TY_B32, TY_B34;
36double 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;
37double 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;
38double TY_w[23];
39
40double TY_sigma( double s, 
41                                 double Z1, double Z2, 
42                             double a, double b, double c1, double c2, double d1, double d2 )
43{
44        return -(a / 2. + b + c1 * exp( -Z1 ) + c2 * exp( -Z2 )) / s + a * pow( s, -3 ) + b * pow( s, -2 ) + 
45        ( c1 + d1 ) * pow( s + Z1, -1 ) + ( c2 + d2 ) * pow( s + Z2, -1 );
46}
47
48double TY_tau( double s, 
49                           double Z1, double Z2, 
50                       double a, double b, double c1, double c2 )
51{
52        return b * pow( s, -2 ) + a * ( pow( s, -3 ) + pow( s, -2 ) ) - pow( s, -1 ) * ( c1 * Z1 * exp( -Z1 ) * 
53                                                                                                                                                                        pow( s + Z1, -1 ) + c2 * Z2 * exp( -Z2 ) * pow( s + Z2, -1 ) );
54} 
55
56double TY_q( double s, 
57                             double Z1, double Z2, 
58                             double a, double b, double c1, double c2, double d1, double d2 )
59{
60        return TY_sigma(s, Z1, Z2, a, b, c1, c2, d1, d2) - exp( -s ) * TY_tau(s, Z1, Z2, a,b, c1, c2);
61}
62
63double TY_g( double s, 
64                     double phi, double Z1, double Z2, 
65                     double a, double b, double c1, double c2, double d1, double d2 )
66{
67        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 ) );
68}
69
70/*
71 ==================================================================================================
72 
73 Structure factor for the potential
74 
75 V(r) = -kB * T * ( K1 * exp[ -Z1 * (r - 1)] / r + K2 * exp[ -Z2 * (r - 1)] / r ) for r > 1
76 V(r) = inf for r <= 1
77 
78 The structure factor is parametrized by (a, b, c1, c2, d1, d2)
79 which depend on (K1, K2, Z1, Z2, phi). 
80 
81 ==================================================================================================
82 */
83
84double TY_hq( double q, double Z, double K, double v )
85{
86        double t1, t2, t3, t4;
87
88        if ( q == 0) 
89        {
90                return (exp(-2*Z)*(v + (v*(-1 + Z) - 2*K*Z)*exp(Z))*(-(v*(1 + Z)) + (v + 2*K*Z*(1 + Z))*exp(Z))*pow(K,-1)*pow(Z,-4))/4.;
91        }
92        else 
93        {
94               
95                t1 = ( 1 - v / ( 2 * K * Z * exp( Z ) ) ) * ( ( 1 - cos( q ) ) / ( q*q ) - 1 / ( Z*Z + q*q ) );
96                t2 = ( v*v * ( q * cos( q ) - Z * sin( q ) ) ) / ( 4 * K * Z*Z * q * ( Z*Z + q*q ) );
97                t3 = ( q * cos( q ) + Z * sin( q ) ) / ( q * ( Z*Z + q*q ) );
98                t4 = v / ( Z * exp( Z ) ) - v*v / ( 4 * K * Z*Z * exp( 2 * Z ) ) - K;
99               
100                return v / Z * t1 - t2 + t3 * t4;
101        }
102}
103
104double TY_pc( double q, 
105                          double Z1, double Z2,double K1, double K2, double phi,
106                          double a, double b, double c1, double c2, double d1, double d2 )
107{       
108        double v1 = 24 * phi * K1 * exp( Z1 ) * TY_g( Z1, phi, Z1, Z2, a, b, c1, c2, d1, d2 );
109        double v2 = 24 * phi * K2 * exp( Z2 ) * TY_g( Z2, phi, Z1, Z2, a, b, c1, c2, d1, d2 );
110       
111        double a0 = a * a;
112        double b0 = -12 * phi *( pow( a + b,2 ) / 2 + a * ( c1 * exp( -Z1 ) + c2 * exp( -Z2 ) ) );
113       
114        double t1, t2, t3,t4;
115       
116        if ( q == 0 ) 
117        {
118                t1 = a0 / 3;
119                t2 = b0 / 4;
120                t3 = a0 * phi / 12;
121        }
122        else 
123        {
124                t1 = a0 * ( sin( q ) - q * cos( q ) ) / pow( q, 3 );
125                t2 = b0 * ( 2 * q * sin( q ) - ( q * q - 2 ) * cos( q ) - 2 ) / pow( q, 4 );
126                t3 = a0 * phi * ( ( q*q - 6 ) * 4 * q * sin( q ) - ( pow( q, 4 ) - 12 * q*q + 24) * cos( q ) + 24 ) / ( 2 * pow( q, 6 ) );
127        }
128        t4 = TY_hq( q, Z1, K1, v1 ) + TY_hq( q, Z2, K2, v2 );
129        return -24 * phi * ( t1 + t2 + t3 + t4 );
130}
131
132double SqTwoYukawa( double q, 
133                                    double Z1, double Z2, double K1, double K2, double phi,
134                                    double a, double b, double c1, double c2, double d1, double d2 )
135{
136        if ( Z1 == Z2 ) 
137        {
138                // one-yukawa potential
139                return 0;
140        } 
141        else 
142        {
143                // two-yukawa potential
144                return 1. / ( 1. - TY_pc( q, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) );
145        }
146}
147
148/*
149==================================================================================================
150
151 Non-linear eqaution system that determines the parameter for structure factor
152 
153==================================================================================================
154*/
155
156double TY_LinearEquation_1( double Z1, double Z2, double K1, double K2, double phi, 
157                                                    double a, double b, double c1, double c2, double d1, double d2 )
158{
159        return b - 12 * phi * ( -a / 8. - b / 6. + d1 * pow( Z1, -2 ) + c1 * ( pow( Z1, -2 )  - exp( -Z1 ) * 
160                    ( 0.5 + ( 1 + Z1 ) * pow( Z1, -2 ) ) ) + d2 * pow( Z2, -2 ) + c2 * ( pow( Z2, -2 ) - exp( -Z2 )
161                        * ( 0.5 + ( 1 + Z2 ) * pow( Z2, -2 ) ) ) );
162}
163
164double TY_LinearEquation_2( double Z1, double Z2, double K1, double K2, double phi, 
165                                                        double a, double b, double c1, double c2, double d1, double d2 )
166{
167        return 1 - a - 12 * phi * ( -a / 3. - b / 2. + d1 * pow( Z1, -1 ) + c1 * ( pow( Z1, -1 ) - ( 1 + Z1 ) * 
168                   exp( -Z1 ) * pow( Z1, -1 ) ) + d2 * pow( Z2, -1 ) + c2 * ( pow( Z2, -1 ) - ( 1 + Z2 ) * exp( -Z2 ) * pow( Z2, -1 ) ) );
169}       
170
171double TY_LinearEquation_3( double Z1, double Z2, double K1, double K2, double phi, 
172                                                    double a, double b, double c1, double c2, double d1, double d2 )
173{                                                               
174        return K1 * exp( Z1 ) - d1 * Z1 * ( 1 - 12 * phi * TY_q( Z1, Z1, Z2, a, b, c1, c2, d1, d2 ) ); 
175}
176
177double TY_LinearEquation_4( double Z1, double Z2, double K1, double K2, double phi, 
178                                                    double a, double b, double c1, double c2, double d1, double d2 )
179{
180        return K2 * exp( Z2 ) - d2 * Z2 * ( 1 - 12 * phi * TY_q( Z2, Z1, Z2, a, b, c1, c2, d1, d2 ) );
181}
182
183double TY_NonlinearEquation_1( double Z1, double Z2, double K1, double K2, double phi, 
184                                                       double a, double b, double c1, double c2, double d1, double d2 )
185{
186        return c1 + d1 - 12 * phi * ( ( c1 + d1 ) * TY_sigma( Z1, Z1, Z2, a, b, c1, c2, d1, d2 ) 
187                                                                 - c1 * TY_tau( Z1, Z1, Z2, a, b, c1, c2 ) * exp( -Z1 ) );
188}
189
190double TY_NonlinearEquation_2( double Z1, double Z2, double K1, double K2, double phi, 
191                                                       double a, double b, double c1, double c2, double d1, double d2 )
192{
193        return c2 + d2 - 12 * phi * ( ( c2 + d2 ) * TY_sigma( Z2, Z1, Z2, a, b, c1, c2, d1, d2 ) 
194                                                                 - c2 * TY_tau( Z2, Z1, Z2, a, b, c1, c2 ) * exp( -Z2 ) );
195}
196
197// Check the computed solutions satisfy the system of equations
198int TY_CheckSolution( double Z1, double Z2, double K1, double K2, double phi, 
199                                      double a, double b, double c1, double c2, double d1, double d2 )
200{
201        double eq_1 = chop( TY_LinearEquation_1( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) );
202        double eq_2 = chop( TY_LinearEquation_2( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) );
203        double eq_3 = chop( TY_LinearEquation_3( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) );
204        double eq_4 = chop( TY_LinearEquation_4( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) );
205        double eq_5 = chop( TY_NonlinearEquation_1( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) );
206        double eq_6 = chop( TY_NonlinearEquation_2( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) );
207       
208        // check if all equation are zero
209        return eq_1 == 0 && eq_2 == 0 && eq_3 == 0 && eq_4 == 0 && eq_5 == 0 && eq_6 == 0;
210}
211
212void TY_ReduceNonlinearSystem( double Z1, double Z2, double K1, double K2, double phi, int debug )
213{
214        /* solution of the 4 linear equations depending on d1 and d2, the solution is polynomial
215         in d1, d2. We represend the solution as determiants obtained by Cramer's rule
216         which can be expressed by their coefficient matrices
217         */
218        char buf[256];
219       
220        double m11 = (3*phi)/2.;
221        double m13 = 6*phi*exp(-Z1)*(2 + Z1*(2 + Z1) - 2*exp(Z1))*pow(Z1,-2);
222        double m14 = 6*phi*exp(-Z2)*(2 + Z2*(2 + Z2) - 2*exp(Z2))*pow(Z2,-2);
223        double m23 = -12*phi*exp(-Z1)*(-1 - Z1 + exp(Z1))*pow(Z1,-1);
224        double m24 = -12*phi*exp(-Z2)*(-1 - Z2 + exp(Z2))*pow(Z2,-1);
225        double m31 = -6*phi*exp(-Z1)*pow(Z1,-2)*(2*(1 + Z1) + exp(Z1)*(-2 + pow(Z1,2)));
226        double m32 = -12*phi*(-1 + Z1 + exp(-Z1))*pow(Z1,-1);
227        double m33 = 6*phi*exp(-2*Z1)*pow(-1 + exp(Z1),2);
228        double m34 = 12*phi*exp(-Z1 - Z2)*(Z2 - (Z1 + Z2)*exp(Z1) + Z1*exp(Z1 + Z2))*pow(Z1 + Z2,-1);
229        double m41 = -6*phi*exp(-Z2)*pow(Z2,-2)*(2*(1 + Z2) + exp(Z2)*(-2 + pow(Z2,2)));
230        double m42 = -12*phi*(-1 + Z2 + exp(-Z2))*pow(Z2,-1);
231        double m43 = 12*phi*exp(-Z1 - Z2)*(Z1 - (Z1 + Z2 - Z2*exp(Z1))*exp(Z2))*pow(Z1 + Z2,-1);       
232        double m44 = 6*phi*exp(-2*Z2)*pow(-1 + exp(Z2),2);
233       
234        /* determinant of the linear system expressed as coefficient matrix in d1, d2 */
235       
236        TY_q22 = m14*(-(m33*m42) + m23*(m32*m41 - m31*m42) + m32*m43 + (4*m11*(-3*m33*m41 + 2*m33*m42 + 3*m31*m43 - 2*m32*m43))/3.) + 
237                  m13*(m34*m42 + m24*(-(m32*m41) + m31*m42) - m32*m44 + (4*m11*(3*m34*m41 - 2*m34*m42 - 3*m31*m44 + 2*m32*m44))/3.) + (3*m24*
238                  (m33*(3*m41 + 4*m11*m41 - 3*m11*m42) + (-3*m31 - 4*m11*m31 + 3*m11*m32)*m43) + 3*m23*(-3*m34*m41 - 4*m11*m34*m41 + 
239                  3*m11*m34*m42 + 3*m31*m44 + 4*m11*m31*m44 - 3*m11*m32*m44) - (m34*m43 - m33*m44)*pow(3 - 2*m11,2))/9.;
240/*     
241        if( debug )
242        {
243                sprintf(buf,"\rDet = \r" );
244                XOPNotice(buf);
245                sprintf(buf, "%f\t%f\r%f\t%f\r", 0., 0., 0., TY_q22 );
246                XOPNotice(buf);
247        }
248*/
249       
250        /* Matrix representation of the determinant of the of the system where row refering to
251         the variable a is replaced by solution vector */
252       
253        TY_qa12 = (K1*(3*m14*(m23*m42 - 4*m11*m43) - 3*m13*(m24*m42 - 4*m11*m44) + (3 + 4*m11)*(m24*m43 - m23*m44))*exp(Z1))/3.;
254       
255        TY_qa21 = -(K2*(3*m14*(m23*m32 - 4*m11*m33) - 3*m13*(m24*m32 - 4*m11*m34) + (3 + 4*m11)*(m24*m33 - m23*m34))*exp(Z2))/3.;
256       
257        TY_qa22 = m14*(-(m23*m42*Z1) + 4*m11*m43*Z1 - m33*(m42 + 4*m11*Z2) + m32*(m43 + m23*Z2)) + 
258                   (3*m13*(m24*m42*Z1 - 4*m11*m44*Z1 + m34*(m42 + 4*m11*Z2) - m32*(m44 + m24*Z2)) + 
259                   (3 + 4*m11)*(-(m24*m43*Z1) + m23*m44*Z1 - m34*(m43 + m23*Z2) + m33*(m44 + m24*Z2)))/3.;     
260       
261        TY_qa23 = 2*phi*pow(Z2,-2)*(m24*(2*(-3*m13*m42 + 3*m43 + 4*m11*m43)*Z1*pow(Z2,2) - m33*(Z1 + Z2)*(6*m42 + (3 + 4*m11)*pow(Z2,2)) + 
262                   3*m32*(Z1 + Z2)*(2*m43 + m13*pow(Z2,2))) + 
263                   m23*(2*(3*m14*m42 - 3*m44 - 4*m11*m44)*Z1*pow(Z2,2) + m34*(Z1 + Z2)*(6*m42 + (3 + 4*m11)*pow(Z2,2)) - 
264                   3*m32*(Z1 + Z2)*(2*m44 + m14*pow(Z2,2))) + 
265                   2*(3*(m14*m33*m42 - m13*m34*m42 - m14*m32*m43 + m34*m43 + m13*m32*m44 - m33*m44)*Z2*(Z1 + Z2) + 
266                   2*m11*(6*(-(m14*m43) + m13*m44)*Z1*pow(Z2,2) + m34*(Z1 + Z2)*(2*m43*(-3 + Z2) - 3*m13*pow(Z2,2)) + 
267                  m33*(Z1 + Z2)*(6*m44 - 2*m44*Z2 + 3*m14*pow(Z2,2)))))*pow(Z1 + Z2,-1);       
268       
269        TY_qa32 = 2*phi*pow(Z1,-2)*(m24*((-3*m13*m42 + (3 + 4*m11)*m43)*(Z1 + Z2)*pow(Z1,2) - 
270                   2*m33*(3*m42*(Z1 + Z2) + (3 + 4*m11)*Z2*pow(Z1,2)) + 6*m32*(m43*(Z1 + Z2) + m13*Z2*pow(Z1,2))) + 
271                   m23*((3*m14*m42 - (3 + 4*m11)*m44)*(Z1 + Z2)*pow(Z1,2) + m34*(6*m42*(Z1 + Z2) + 2*(3 + 4*m11)*Z2*pow(Z1,2)) - 
272                   6*m32*(m44*(Z1 + Z2) + m14*Z2*pow(Z1,2))) + 
273                   2*(3*(m14*m33*m42 - m13*m34*m42 - m14*m32*m43 + m34*m43 + m13*m32*m44 - m33*m44)*Z1*(Z1 + Z2) + 
274                   2*m11*(-3*(m14*m43 - m13*m44)*(Z1 + Z2)*pow(Z1,2) + 2*m34*(m43*(-3 + Z1)*(Z1 + Z2) - 3*m13*Z2*pow(Z1,2)) + 
275                  m33*(-2*m44*(-3 + Z1)*(Z1 + Z2) + 6*m14*Z2*pow(Z1,2)))))*pow(Z1 + Z2,-1);
276/*             
277        if( debug )
278        {
279                sprintf(buf,"\rDet_a = \r" );
280                XOPNotice(buf);
281                sprintf(buf, "%f\t%f\t%f\r%f\t%f\t%f\r%f\t%f\t%f\r",
282                           0., TY_qa12, 0.,
283                           TY_qa21, TY_qa22, TY_qa23,
284                           0.,    TY_qa32, 0. );
285                XOPNotice(buf);
286        }
287*/
288       
289        /* Matrix representation of the determinant of the of the system where row refering to
290         the variable b is replaced by solution vector */
291       
292        TY_qb12 = (K1*(-3*m11*m24*m43 + m14*(-3*m23*m41 + (-3 + 8*m11)*m43) + 3*m11*m23*m44 + m13*(3*m24*m41 + 3*m44 - 8*m11*m44))*exp(Z1))/3.;
293       
294        TY_qb21 = (K2*(-3*m13*m24*m31 + 3*m11*m24*m33 + m14*(3*m23*m31 + (3 - 8*m11)*m33) - 3*m13*m34 + 8*m11*m13*m34 - 3*m11*m23*m34)*
295                        exp(Z2))/3.;
296       
297        TY_qb22 = m13*(m31*m44 - m24*m41*Z1 - m44*Z1 + (8*m11*m44*Z1)/3. + m24*m31*Z2 + m34*(-m41 + Z2 - (8*m11*Z2)/3.)) + 
298                   m14*(m23*m41*Z1 + m43*Z1 - (8*m11*m43*Z1)/3. + m33*(m41 - Z2 + (8*m11*Z2)/3.) - m31*(m43 + m23*Z2)) + 
299                   m11*(m24*m43*Z1 - m23*m44*Z1 + m34*(m43 + m23*Z2) - m33*(m44 + m24*Z2));     
300       
301        TY_qb23 = 2*phi*(3*m14*m23*m31 - 3*m13*m24*m31 + 3*m14*m33 - 8*m11*m14*m33 + 3*m11*m24*m33 - 3*m13*m34 + 8*m11*m13*m34 - 
302                   3*m11*m23*m34 + 2*(3*m24*(m33*m41 - m31*m43) + m23*(-3*m34*m41 + 3*m31*m44) + (-3 + 8*m11)*(m34*m43 - m33*m44))*
303                   pow(Z2,-2) + 6*(-(m14*m33*m41) + m13*m34*m41 + m14*m31*m43 - m11*m34*m43 - m13*m31*m44 + m11*m33*m44)*pow(Z2,-1) + 
304                   2*(-3*m11*m24*m43 + m14*(-3*m23*m41 + (-3 + 8*m11)*m43) + 3*m11*m23*m44 + m13*(3*m24*m41 + 3*m44 - 8*m11*m44))*Z1*
305                   pow(Z1 + Z2,-1));
306       
307        TY_qb32 = 2*phi*pow(Z1,-2)*(6*(-(m34*(m23*m41 + m43)) + m24*(m33*m41 - m31*m43) + (m23*m31 + m33)*m44) + 
308                   6*(-(m14*m33*m41) + m13*m34*m41 + m14*m31*m43 - m13*m31*m44)*Z1 + 
309                   3*(m14*(2*m23*m31 + 2*m33 - m23*m41 - m43) + m13*(-2*m34 + m24*(-2*m31 + m41) + m44))*pow(Z1,2) + 
310                   (m11*Z2*(16*m34*m43 - 16*m33*m44 - 6*m34*m43*Z1 + 6*m33*m44*Z1 + 
311                   (6*m24*m33 - 3*m24*m43 + 8*m14*(-2*m33 + m43) + (8*m13 - 3*m23)*(2*m34 - m44))*pow(Z1,2)) + 
312                   m11*Z1*(2*m34*m43*(8 - 3*Z1) + 2*m33*m44*(-8 + 3*Z1) + (8*m14*m43 - 3*m24*m43 - 8*m13*m44 + 3*m23*m44)*pow(Z1,2)))*
313                   pow(Z1 + Z2,-1) + 6*(-(m14*(m23*m31 + m33)) + m13*(m24*m31 + m34))*pow(Z1,3)*pow(Z1 + Z2,-1));
314/*             
315        if( debug )
316        {
317                sprintf(buf,"\rDet_b = \r" );
318                XOPNotice(buf);
319                sprintf(buf, "%f\t%f\t%f\r%f\t%f\t%f\r%f\t%f\t%f\r",
320                           0., TY_qb12, 0.,
321                           TY_qb21, TY_qb22, TY_qb23,
322                           0., TY_qb32, 0. );
323                XOPNotice(buf);
324        }
325*/
326       
327        /* Matrix representation of the determinant of the of the system where row refering to
328         the variable c1 is replaced by solution vector */
329       
330        TY_qc112 = -(K1*exp(Z1)*(9*m24*m41 - 9*m14*m42 + 3*m11*(-12*m14*m41 + 4*m24*m41 + 8*m14*m42 - 3*m24*m42) + m44*pow(3 - 2*m11,2)))/9.;
331       
332        TY_qc121 = (K2*exp(Z2)*(9*m24*m31 - 9*m14*m32 + 3*m11*(-12*m14*m31 + 4*m24*m31 + 8*m14*m32 - 3*m24*m32) + m34*pow(3 - 2*m11,2)))/9.;
333       
334        TY_qc122 = m14*(-4*m11*m41*Z1 - m42*Z1 + (8*m11*m42*Z1)/3. + m32*(-m41 + Z2 - (8*m11*Z2)/3.) + m31*(m42 + 4*m11*Z2)) + 
335                        (3*m34*((3 + 4*m11)*m41 - 3*m11*m42) + 9*m11*m32*m44 + 9*m24*m41*Z1 + 12*m11*m24*m41*Z1 - 9*m11*m24*m42*Z1 + 9*m44*Z1 - 
336                        12*m11*m44*Z1 + 9*m11*m24*m32*Z2 - 3*(3 + 4*m11)*m31*(m44 + m24*Z2) - m34*Z2*pow(3 - 2*m11,2) + 4*m44*Z1*pow(m11,2))/9.;
337       
338        TY_qc123 = (2*phi*pow(Z2,-2)*(9*(m34*(Z1 + Z2)*(2*m42 + Z2*(-2*m41 + Z2)) - m32*(Z1 + Z2)*(2*m44 + m14*Z2*(-2*m41 + Z2)) - 
339                        2*(m14*m42 - m44)*Z2*(-(Z1*Z2) + m31*(Z1 + Z2))) + 4*(-2*m44*Z1 + m34*(Z1 + Z2))*pow(m11,2)*pow(Z2,2) - 
340                        3*m24*(2*(3*m41 + 4*m11*m41 - 3*m11*m42)*Z1*pow(Z2,2) + 3*m32*(Z1 + Z2)*(2*m41 + m11*pow(Z2,2)) - 
341                        m31*(Z1 + Z2)*(6*m42 + (3 + 4*m11)*pow(Z2,2))) - 
342                        6*m11*(-8*m32*m44*Z1 + m32*m44*(-8 + 3*Z1)*Z2 + (3*m32*m44 - 4*(m14*(m32 + 3*m41 - 2*m42) + m44)*Z1)*pow(Z2,2) + 
343                        m34*(Z1 + Z2)*(8*m42 + 4*m41*(-3 + Z2) - 3*m42*Z2 + 2*pow(Z2,2)) + 2*m31*(Z1 + Z2)*(6*m44 - 2*m44*Z2 + 3*m14*pow(Z2,2)) - 
344                        4*m14*m32*pow(Z2,3)))*pow(Z1 + Z2,-1))/3.;
345       
346        TY_qc132 = (-2*phi*pow(Z1,-2)*(9*((m14*m42 - m44)*(2*m31 - Z1)*Z1*(Z1 + Z2) - 2*m34*(m42*(Z1 + Z2) - Z1*(-(Z1*Z2) + m41*(Z1 + Z2))) + 
347                        2*m32*(m44*(Z1 + Z2) - m14*Z1*(-(Z1*Z2) + m41*(Z1 + Z2)))) + 4*(-2*m34*Z2 + m44*(Z1 + Z2))*pow(m11,2)*pow(Z1,2) + 
348                        3*m24*(((3 + 4*m11)*m41 - 3*m11*m42)*(Z1 + Z2)*pow(Z1,2) + 6*m32*(m41*(Z1 + Z2) + m11*Z2*pow(Z1,2)) - 
349                        2*m31*(3*m42*(Z1 + Z2) + (3 + 4*m11)*Z2*pow(Z1,2))) + 
350                        6*m11*(Z1*(-8*m32*m44 + m34*(m42*(8 - 3*Z1) + 4*m41*(-3 + Z1)) - 4*m31*m44*(-3 + Z1) + 3*m32*m44*Z1 - 
351                        2*(3*m14*m41 - 2*m14*m42 + m44)*pow(Z1,2)) + 
352                        Z2*(4*(3*m31 - 2*m32)*m44 + Z1*(-4*m31*m44 + 3*m32*m44 - 2*(m14*(-6*m31 + 4*m32 + 3*m41 - 2*m42) + m44)*Z1) + 
353                        m34*(m42*(8 - 3*Z1) + 4*m41*(-3 + Z1) + 4*pow(Z1,2)))))*pow(Z1 + Z2,-1))/3.;
354/*             
355        if( debug )
356        {
357                sprintf(buf,"\rDet_c1 = \r" );
358                XOPNotice(buf);
359                sprintf(buf, "%f\t%f\t%f\r%f\t%f\t%f\r%f\t%f\t%f\r",
360                            0., TY_qc112, 0.,
361                           TY_qc121, TY_qc122, TY_qc123,
362                           0., TY_qc132, 0. );
363                XOPNotice(buf);
364        }
365 */
366        /* Matrix representation of the determinant of the of the system where row refering to
367         the variable c1 is replaced by solution vector */
368       
369        TY_qc212 = (K1*exp(Z1)*(9*m23*m41 - 9*m13*m42 + 3*m11*(-12*m13*m41 + 4*m23*m41 + 8*m13*m42 - 3*m23*m42) + m43*pow(3 - 2*m11,2)))/9.;
370       
371        TY_qc221 = -(K2*exp(Z2)*(9*m23*m31 - 9*m13*m32 + 3*m11*(-12*m13*m31 + 4*m23*m31 + 8*m13*m32 - 3*m23*m32) + m33*pow(3 - 2*m11,2)))/9.;
372       
373        TY_qc222 = m13*(4*m11*m41*Z1 + m42*Z1 - (8*m11*m42*Z1)/3. + m32*(m41 - Z2 + (8*m11*Z2)/3.) - m31*(m42 + 4*m11*Z2)) + 
374                        (9*m31*m43 - 9*(m23*m41 + m43)*Z1 + 9*m23*m31*Z2 + 3*m11*
375                        ((-4*m23*m41 + 3*m23*m42 + 4*m43)*Z1 + 4*m31*(m43 + m23*Z2) - 3*m32*(m43 + m23*Z2)) + 
376                        m33*(-3*(3 + 4*m11)*m41 + 9*m11*m42 + Z2*pow(3 - 2*m11,2)) - 4*m43*Z1*pow(m11,2))/9.;
377       
378        TY_qc223 = (2*phi*pow(Z2,-2)*(9*(-(m33*(Z1 + Z2)*(2*m42 + Z2*(-2*m41 + Z2))) + m32*(Z1 + Z2)*(2*m43 + m13*Z2*(-2*m41 + Z2)) + 
379                        2*(m13*m42 - m43)*Z2*(-(Z1*Z2) + m31*(Z1 + Z2))) - 4*(-2*m43*Z1 + m33*(Z1 + Z2))*pow(m11,2)*pow(Z2,2) + 
380                        3*m23*(2*(3*m41 + 4*m11*m41 - 3*m11*m42)*Z1*pow(Z2,2) + 3*m32*(Z1 + Z2)*(2*m41 + m11*pow(Z2,2)) - 
381                        m31*(Z1 + Z2)*(6*m42 + (3 + 4*m11)*pow(Z2,2))) + 
382                        6*m11*(-8*m32*m43*Z1 + m32*m43*(-8 + 3*Z1)*Z2 + (3*m32*m43 - 4*(m13*(m32 + 3*m41 - 2*m42) + m43)*Z1)*pow(Z2,2) + 
383                        m33*(Z1 + Z2)*(8*m42 + 4*m41*(-3 + Z2) - 3*m42*Z2 + 2*pow(Z2,2)) + 2*m31*(Z1 + Z2)*(6*m43 - 2*m43*Z2 + 3*m13*pow(Z2,2)) - 
384                        4*m13*m32*pow(Z2,3)))*pow(Z1 + Z2,-1))/3.;
385       
386        TY_qc232 = (2*phi*pow(Z1,-2)*(9*((m13*m42 - m43)*(2*m31 - Z1)*Z1*(Z1 + Z2) - 2*m33*(m42*(Z1 + Z2) - Z1*(-(Z1*Z2) + m41*(Z1 + Z2))) + 
387                        2*m32*(m43*(Z1 + Z2) - m13*Z1*(-(Z1*Z2) + m41*(Z1 + Z2)))) + 4*(-2*m33*Z2 + m43*(Z1 + Z2))*pow(m11,2)*pow(Z1,2) + 
388                        3*m23*(((3 + 4*m11)*m41 - 3*m11*m42)*(Z1 + Z2)*pow(Z1,2) + 6*m32*(m41*(Z1 + Z2) + m11*Z2*pow(Z1,2)) - 
389                        2*m31*(3*m42*(Z1 + Z2) + (3 + 4*m11)*Z2*pow(Z1,2))) + 
390                        6*m11*(Z1*(-8*m32*m43 + m33*(m42*(8 - 3*Z1) + 4*m41*(-3 + Z1)) - 4*m31*m43*(-3 + Z1) + 3*m32*m43*Z1 - 
391                        2*(3*m13*m41 - 2*m13*m42 + m43)*pow(Z1,2)) + 
392                        Z2*(4*(3*m31 - 2*m32)*m43 + Z1*(-4*m31*m43 + 3*m32*m43 - 2*(m13*(-6*m31 + 4*m32 + 3*m41 - 2*m42) + m43)*Z1) + 
393                        m33*(m42*(8 - 3*Z1) + 4*m41*(-3 + Z1) + 4*pow(Z1,2)))))*pow(Z1 + Z2,-1))/3.;
394/*             
395        if( debug )
396        {
397                sprintf(buf,"\rDet_c2 = \r" );
398                XOPNotice(buf);
399                sprintf(buf, "%f\t%f\t%f\r%f\t%f\t%f\r%f\t%f\t%f\r",
400                           0., TY_qc212, 0.,
401                           TY_qc221, TY_qc222, TY_qc223,
402                           0., TY_qc232, 0. );
403                XOPNotice(buf);
404        }       
405*/
406       
407        /* coefficient matrices of nonlinear equation 1 */
408       
409        TY_A12 = 6*phi*TY_qc112*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*TY_qc212*exp(Z1)*(-Z2 + (Z1 + Z2)*exp(Z1))*pow(Z1,2) + 
410                  exp(Z2)*(exp(2*Z1)*(Z1*(2*TY_qb12*(-1 + Z1)*(Z1 + Z2) - Z1*(2*TY_qc212*Z1 + TY_qc112*(Z1 + Z2))) + TY_qa12*(Z1 + Z2)*(-2 + pow(Z1,2))) - 
411                  TY_qc112*(Z1 + Z2)*pow(Z1,2) + 2*(Z1 + Z2)*exp(Z1)*(TY_qa12 + (TY_qa12 + TY_qb12)*Z1 + TY_qc112*pow(Z1,2))))*pow(Z1 + Z2,-1);
412       
413        TY_A21 = 6*phi*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*(TY_qc121*TY_qc212 + TY_qc112*TY_qc221)*exp(Z1)*(-Z2 + (Z1 + Z2)*exp(Z1))*pow(Z1,2) + 
414                  exp(Z2)*(2*(TY_qa12*TY_qc121 + TY_qa21*TY_qc112*(1 + Z1) + Z1*(TY_qb21*TY_qc112 + TY_qc121*(TY_qa12 + TY_qb12 + 2*TY_qc112*Z1)))*(Z1 + Z2)*exp(Z1) + 
415                  exp(2*Z1)*(2*Z1*(TY_qb21*TY_qc112*(-1 + Z1)*(Z1 + Z2) + TY_qb12*TY_qc121*(-1 + Z1)*(Z1 + Z2) - 
416                  Z1*(TY_qc121*TY_qc212*Z1 + TY_qc112*(TY_qc121 + TY_qc221)*Z1 + TY_qc112*TY_qc121*Z2)) + TY_qa21*TY_qc112*(Z1 + Z2)*(-2 + pow(Z1,2)) + 
417                  TY_qa12*TY_qc121*(Z1 + Z2)*(-2 + pow(Z1,2))) - 2*TY_qc112*TY_qc121*(Z1 + Z2)*pow(Z1,2)))*pow(Z1 + Z2,-1);
418       
419        TY_A22 = exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(12*phi*(TY_qc122*TY_qc212 + TY_qc112*TY_qc222)*exp(Z1)*(-Z2 + (Z1 + Z2)*exp(Z1))*pow(Z1,2) + 
420                  exp(Z2)*(12*phi*(TY_qa12*TY_qc122 + TY_qa22*TY_qc112*(1 + Z1) + Z1*(TY_qb22*TY_qc112 + TY_qc122*(TY_qa12 + TY_qb12 + 2*TY_qc112*Z1)))*(Z1 + Z2)*exp(Z1) - 
421                  2*phi*TY_qc112*TY_qc122*(Z1 + Z2)*pow(Z1,2) + exp(2*Z1)*
422                  (6*phi*(2*Z1*(TY_qb22*TY_qc112*(-1 + Z1)*(Z1 + Z2) + TY_qb12*TY_qc122*(-1 + Z1)*(Z1 + Z2) - 
423                  Z1*(TY_qc122*TY_qc212*Z1 + TY_qc112*(TY_qc122 + TY_qc222)*Z1 + TY_qc112*TY_qc122*Z2)) + TY_qa22*TY_qc112*(Z1 + Z2)*(-2 + pow(Z1,2)) + 
424                  TY_qa12*TY_qc122*(Z1 + Z2)*(-2 + pow(Z1,2))) + TY_q22*TY_qc112*(Z1 + Z2)*pow(Z1,3))))*pow(Z1 + Z2,-1);
425       
426        TY_A23 = 6*phi*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*(TY_qc123*TY_qc212 + TY_qc112*TY_qc223)*exp(Z1)*(-Z2 + (Z1 + Z2)*exp(Z1))*pow(Z1,2) + 
427                  exp(Z2)*(2*(TY_qa12*TY_qc123 + TY_qa23*TY_qc112*(1 + Z1) + Z1*(TY_qb23*TY_qc112 + TY_qc123*(TY_qa12 + TY_qb12 + 2*TY_qc112*Z1)))*(Z1 + Z2)*exp(Z1) + 
428                  exp(2*Z1)*(2*Z1*(TY_qb23*TY_qc112*(-1 + Z1)*(Z1 + Z2) + TY_qb12*TY_qc123*(-1 + Z1)*(Z1 + Z2) - 
429                  Z1*((TY_q22*TY_qc112 + TY_qc123*(TY_qc112 + TY_qc212) + TY_qc112*TY_qc223)*Z1 + TY_qc112*TY_qc123*Z2)) + TY_qa23*TY_qc112*(Z1 + Z2)*(-2 + pow(Z1,2)) + 
430                  TY_qa12*TY_qc123*(Z1 + Z2)*(-2 + pow(Z1,2))) - 2*TY_qc112*TY_qc123*(Z1 + Z2)*pow(Z1,2)))*pow(Z1 + Z2,-1);
431       
432        TY_A32 = exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(12*phi*(TY_qa23*TY_qc121 + TY_qa22*TY_qc122 + TY_qa21*TY_qc123 + TY_qa12*TY_qc132 + TY_qa32*TY_qc112*(1 + Z1) + 
433                  Z1*(TY_qb32*TY_qc112 + (TY_qa23 + TY_qb23)*TY_qc121 + (TY_qa21 + TY_qb21)*TY_qc123 + (TY_qa12 + TY_qb12)*TY_qc132 + TY_q22*TY_qc112*Z1 + 
434                  2*(TY_qc121*TY_qc123 + TY_qc112*TY_qc132)*Z1 + TY_qc122*(TY_qa22 + TY_qb22 + TY_qc122*Z1)))*(Z1 + Z2)*exp(Z1 + Z2) - 
435                  12*phi*(TY_qc132*TY_qc212 + TY_qc123*TY_qc221 + TY_qc122*TY_qc222 + TY_qc121*TY_qc223 + TY_qc112*TY_qc232)*Z2*exp(Z1)*pow(Z1,2) + 
436                  12*phi*((TY_q22 + TY_qc132)*TY_qc212 + TY_qc123*TY_qc221 + TY_qc122*TY_qc222 + TY_qc121*TY_qc223 + TY_qc112*TY_qc232)*(Z1 + Z2)*exp(2*Z1)*pow(Z1,2) - 
437                  6*phi*(Z1 + Z2)*exp(Z2)*(2*TY_qc121*TY_qc123 + 2*TY_qc112*TY_qc132 + pow(TY_qc122,2))*pow(Z1,2) + 
438                  exp(2*Z1 + Z2)*(TY_q22*(6*phi*(2*Z1*(TY_qb12*(-1 + Z1)*(Z1 + Z2) - Z1*((TY_qc112 + TY_qc121 + TY_qc212)*Z1 + TY_qc112*Z2)) + 
439                  TY_qa12*(Z1 + Z2)*(-2 + pow(Z1,2))) + TY_qc122*(Z1 + Z2)*pow(Z1,3)) + 
440                  6*phi*(-2*TY_qa22*TY_qc122*Z1 - 2*TY_qa21*TY_qc123*Z1 - 2*TY_qa12*TY_qc132*Z1 + TY_qa32*TY_qc112*(Z1 + Z2)*(-2 + pow(Z1,2)) + 
441                  TY_qa23*TY_qc121*(Z1 + Z2)*(-2 + pow(Z1,2)) - 2*TY_qb32*TY_qc112*pow(Z1,2) - 2*TY_qb23*TY_qc121*pow(Z1,2) - 2*TY_qb22*TY_qc122*pow(Z1,2) - 
442                  2*TY_qb21*TY_qc123*pow(Z1,2) - 2*TY_qb12*TY_qc132*pow(Z1,2) + 
443                  Z2*(-2*(TY_qa22*TY_qc122 + TY_qa21*TY_qc123 + TY_qa12*TY_qc132) - 2*(TY_qb32*TY_qc112 + TY_qb23*TY_qc121 + TY_qb22*TY_qc122 + TY_qb21*TY_qc123 + TY_qb12*TY_qc132)*Z1 + 
444                  (2*TY_qb32*TY_qc112 + 2*TY_qb23*TY_qc121 + (TY_qa22 + 2*TY_qb22 - TY_qc122)*TY_qc122 + (TY_qa21 + 2*TY_qb21 - 2*TY_qc121)*TY_qc123 + 
445                  (TY_qa12 + 2*TY_qb12 - 2*TY_qc112)*TY_qc132)*pow(Z1,2)) + 2*TY_qb32*TY_qc112*pow(Z1,3) + 2*TY_qb23*TY_qc121*pow(Z1,3) + TY_qa22*TY_qc122*pow(Z1,3) + 
446                  2*TY_qb22*TY_qc122*pow(Z1,3) + TY_qa21*TY_qc123*pow(Z1,3) + 2*TY_qb21*TY_qc123*pow(Z1,3) - 2*TY_qc121*TY_qc123*pow(Z1,3) + TY_qa12*TY_qc132*pow(Z1,3) + 
447                  2*TY_qb12*TY_qc132*pow(Z1,3) - 2*TY_qc112*TY_qc132*pow(Z1,3) - 2*TY_qc132*TY_qc212*pow(Z1,3) - 2*TY_qc123*TY_qc221*pow(Z1,3) - 
448                  2*TY_qc122*TY_qc222*pow(Z1,3) - 2*TY_qc121*TY_qc223*pow(Z1,3) - 2*TY_qc112*TY_qc232*pow(Z1,3) - pow(TY_qc122,2)*pow(Z1,3))))*pow(Z1 + Z2,-1);
449       
450        TY_A41 = 6*phi*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*exp(Z1)*
451                  ((-(TY_qc132*TY_qc221) - TY_qc121*TY_qc232)*Z2 + ((TY_q22 + TY_qc132)*TY_qc221 + TY_qc121*TY_qc232)*(Z1 + Z2)*exp(Z1))*pow(Z1,2) + 
452                  exp(Z2)*(2*(TY_qa21*TY_qc132 + TY_qa32*TY_qc121*(1 + Z1) + Z1*(TY_qb32*TY_qc121 + (TY_qa21 + TY_qb21)*TY_qc132 + TY_qc121*(TY_q22 + 2*TY_qc132)*Z1))*(Z1 + Z2)*
453                  exp(Z1) - 2*TY_qc121*TY_qc132*(Z1 + Z2)*pow(Z1,2) + 
454                  exp(2*Z1)*(Z2*(-2*(TY_qa32*TY_qc121 + TY_qa21*(TY_q22 + TY_qc132)) - 2*(TY_qb32*TY_qc121 + TY_qb21*(TY_q22 + TY_qc132))*Z1 + 
455                  (TY_q22*(TY_qa21 + 2*TY_qb21 - 2*TY_qc121) + TY_qc121*(TY_qa32 + 2*TY_qb32 - 2*TY_qc132) + (TY_qa21 + 2*TY_qb21)*TY_qc132)*pow(Z1,2)) + 
456                  Z1*(-2*(TY_qa32*TY_qc121 + TY_qa21*(TY_q22 + TY_qc132)) - 2*(TY_qb32*TY_qc121 + TY_qb21*(TY_q22 + TY_qc132))*Z1 + 
457                  (TY_qa32*TY_qc121 + 2*TY_qb32*TY_qc121 + TY_qa21*TY_qc132 + 2*TY_qb21*TY_qc132 - 2*TY_qc121*TY_qc132 + TY_q22*(TY_qa21 + 2*TY_qb21 - 2*TY_qc121 - 2*TY_qc221) - 
458                  2*TY_qc132*TY_qc221 - 2*TY_qc121*TY_qc232)*pow(Z1,2)))))*pow(Z1 + Z2,-1);
459       
460        TY_A42 = exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(12*phi*(TY_qa22*TY_qc132 + TY_qa32*TY_qc122*(1 + Z1) + 
461                  Z1*(TY_qb32*TY_qc122 + (TY_qa22 + TY_qb22)*TY_qc132 + TY_qc122*(TY_q22 + 2*TY_qc132)*Z1))*(Z1 + Z2)*exp(Z1 + Z2) - 
462                  12*phi*(TY_qc132*TY_qc222 + TY_qc122*TY_qc232)*Z2*exp(Z1)*pow(Z1,2) + 
463                  12*phi*((TY_q22 + TY_qc132)*TY_qc222 + TY_qc122*TY_qc232)*(Z1 + Z2)*exp(2*Z1)*pow(Z1,2) - 12*phi*TY_qc122*TY_qc132*(Z1 + Z2)*exp(Z2)*pow(Z1,2) + 
464                  exp(2*Z1 + Z2)*(6*phi*(2*Z1*(TY_qb32*TY_qc122*(-1 + Z1)*(Z1 + Z2) + TY_qb22*TY_qc132*(-1 + Z1)*(Z1 + Z2) - 
465                  Z1*(TY_qc132*TY_qc222*Z1 + TY_qc122*(TY_qc132 + TY_qc232)*Z1 + TY_qc122*TY_qc132*Z2)) + TY_qa32*TY_qc122*(Z1 + Z2)*(-2 + pow(Z1,2)) + 
466                  TY_qa22*TY_qc132*(Z1 + Z2)*(-2 + pow(Z1,2))) + (Z1 + Z2)*pow(TY_q22,2)*pow(Z1,3) + 
467                  TY_q22*(6*phi*(2*Z1*(TY_qb22*(-1 + Z1)*(Z1 + Z2) - Z1*((TY_qc122 + TY_qc222)*Z1 + TY_qc122*Z2)) + TY_qa22*(Z1 + Z2)*(-2 + pow(Z1,2))) + 
468                  TY_qc132*(Z1 + Z2)*pow(Z1,3))))*pow(Z1 + Z2,-1);
469       
470        TY_A43 = -6*phi*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*exp(Z1)*
471                  ((TY_qc132*TY_qc223 + TY_qc123*TY_qc232)*Z2 - ((TY_q22 + TY_qc132)*TY_qc223 + TY_qc123*TY_qc232)*(Z1 + Z2)*exp(Z1))*pow(Z1,2) + 
472                  exp(Z2)*(-2*(TY_qa23*TY_qc132 + TY_qa32*TY_qc123*(1 + Z1) + Z1*(TY_qb32*TY_qc123 + (TY_qa23 + TY_qb23)*TY_qc132 + TY_qc123*(TY_q22 + 2*TY_qc132)*Z1))*(Z1 + Z2)*
473                  exp(Z1) + 2*TY_qc123*TY_qc132*(Z1 + Z2)*pow(Z1,2) + 
474                  exp(2*Z1)*(Z2*(2*TY_qa32*TY_qc123 + 2*TY_qa23*(TY_q22 + TY_qc132) + 2*(TY_qb32*TY_qc123 + TY_qb23*(TY_q22 + TY_qc132))*Z1 - 
475                  (TY_q22*(TY_qa23 + 2*TY_qb23 - 2*TY_qc123) + TY_qc123*(TY_qa32 + 2*TY_qb32 - 2*TY_qc132) + (TY_qa23 + 2*TY_qb23)*TY_qc132)*pow(Z1,2)) + 
476                  Z1*(2*TY_qa32*TY_qc123 + 2*TY_qa23*(TY_q22 + TY_qc132) + 2*(TY_qb32*TY_qc123 + TY_qb23*(TY_q22 + TY_qc132))*Z1 + 
477                  (-(TY_qa32*TY_qc123) - (TY_qa23 + 2*TY_qb23)*TY_qc132 + TY_q22*(-TY_qa23 + 2*(-TY_qb23 + TY_qc123 + TY_qc132 + TY_qc223)) + 
478                  2*(-(TY_qb32*TY_qc123) + TY_qc132*(TY_qc123 + TY_qc223) + TY_qc123*TY_qc232) + 2*pow(TY_q22,2))*pow(Z1,2)))))*pow(Z1 + Z2,-1);
479       
480        TY_A52 = -6*phi*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*TY_qc232*exp(Z1)*(TY_qc132*Z2 - (TY_q22 + TY_qc132)*(Z1 + Z2)*exp(Z1))*pow(Z1,2) + 
481                  exp(Z2)*((TY_q22 + TY_qc132)*exp(2*Z1)*(Z1*(-2*TY_qb32*(-1 + Z1)*(Z1 + Z2) + Z1*((TY_q22 + TY_qc132 + 2*TY_qc232)*Z1 + (TY_q22 + TY_qc132)*Z2)) - 
482                  TY_qa32*(Z1 + Z2)*(-2 + pow(Z1,2))) + (Z1 + Z2)*pow(TY_qc132,2)*pow(Z1,2) - 
483                  2*TY_qc132*(Z1 + Z2)*exp(Z1)*(TY_qa32 + (TY_qa32 + TY_qb32)*Z1 + (TY_q22 + TY_qc132)*pow(Z1,2))))*pow(Z1 + Z2,-1);
484       
485        // normalize A
486        /*double norm_A = sqrt(pow(TY_A52,2)+pow(TY_A43,2)+pow(TY_A42, 2)+pow(TY_A41, 2)+pow(TY_A32, 2)+
487                                                 pow(TY_A23,2)+pow(TY_A22,2)+pow(TY_A21, 2)+pow(TY_A12, 2));
488        TY_A12 /= norm_A;
489        TY_A21 /= norm_A;
490        TY_A22 /= norm_A;
491        TY_A23 /= norm_A;
492        TY_A32 /= norm_A;
493        TY_A41 /= norm_A;
494        TY_A42 /= norm_A;
495        TY_A43 /= norm_A;
496        TY_A52 /= norm_A;*/
497/*     
498        if( debug )
499        {
500                sprintf(buf,"\rNonlinear equation 1 = \r" );
501                XOPNotice(buf);
502                sprintf(buf, "%f\t\t%f\t\t%f\r", 0.,   TY_A12, 0.   );
503                XOPNotice(buf);
504                sprintf(buf, "%f\t\t%f\t\t%f\r", TY_A21, TY_A22, TY_A23 );
505                XOPNotice(buf);
506                sprintf(buf, "%f\t\t%f\t\t%f\r",  0.,  TY_A32, 0.   );
507                XOPNotice(buf);
508                sprintf(buf, "%f\t\t%f\t\t%f\r", TY_A41, TY_A42, TY_A43 );
509                XOPNotice(buf);
510                sprintf(buf, "%f\t\t%f\t\t%f\r", 0.,   TY_A52, 0. );           
511                XOPNotice(buf);
512        }
513 */
514       
515        /* coefficient matrices of nonlinear equation 2 */
516       
517        TY_B12 = 6*phi*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-2*TY_qc212*TY_qc221*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 
518                  2*exp(Z2)*((Z1 + Z2)*(TY_qa12*TY_qc221 + TY_qa21*TY_qc212*(1 + Z2) + Z2*(TY_qb21*TY_qc212 + TY_qc221*(TY_qa12 + TY_qb12 + 2*TY_qc212*Z2)))*exp(Z1) + 
519                  (-(TY_qc121*TY_qc212) - TY_qc112*TY_qc221)*Z1*pow(Z2,2)) + 
520                  exp(2*Z2)*(exp(Z1)*(2*Z2*(TY_qb21*TY_qc212*(-1 + Z2)*(Z1 + Z2) + TY_qb12*TY_qc221*(-1 + Z2)*(Z1 + Z2) - 
521                  Z2*(TY_qc212*TY_qc221*Z1 + TY_qc112*TY_qc221*Z2 + TY_qc212*(TY_qc121 + TY_qc221)*Z2)) + TY_qa21*TY_qc212*(Z1 + Z2)*(-2 + pow(Z2,2)) + 
522                  TY_qa12*TY_qc221*(Z1 + Z2)*(-2 + pow(Z2,2))) + 2*(TY_qc121*TY_qc212 + TY_qc112*TY_qc221)*(Z1 + Z2)*pow(Z2,2)))*pow(Z1 + Z2,-1);
523       
524        TY_B14 = 6*phi*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-2*TY_qc212*TY_qc223*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 
525                  2*exp(Z2)*((Z1 + Z2)*(TY_qa12*TY_qc223 + TY_qa23*TY_qc212*(1 + Z2) + Z2*(TY_qb23*TY_qc212 + (TY_qa12 + TY_qb12)*TY_qc223 + TY_qc212*(TY_q22 + 2*TY_qc223)*Z2))*
526                  exp(Z1) + (-(TY_qc123*TY_qc212) - TY_qc112*TY_qc223)*Z1*pow(Z2,2)) + 
527                  exp(2*Z2)*(2*(TY_qc123*TY_qc212 + TY_qc112*(TY_q22 + TY_qc223))*(Z1 + Z2)*pow(Z2,2) + 
528                  exp(Z1)*(-2*(TY_qa23*TY_qc212 + TY_qa12*(TY_q22 + TY_qc223))*Z1 - 
529                  2*(TY_qa23*TY_qc212 + TY_qa12*(TY_q22 + TY_qc223) + (TY_qb23*TY_qc212 + TY_qb12*(TY_q22 + TY_qc223))*Z1)*Z2 + 
530                  (-2*(TY_qb23*TY_qc212 + TY_qb12*(TY_q22 + TY_qc223)) + (TY_q22*(TY_qa12 + 2*TY_qb12 - 2*TY_qc212) + TY_qc212*(TY_qa23 + 2*TY_qb23 - 2*TY_qc223) + 
531                  (TY_qa12 + 2*TY_qb12)*TY_qc223)*Z1)*pow(Z2,2) + 
532                  (TY_q22*(TY_qa12 + 2*TY_qb12 - 2*TY_qc112 - 2*TY_qc212) + TY_qc212*(TY_qa23 + 2*TY_qb23 - 2*TY_qc123 - 2*TY_qc223) + (TY_qa12 + 2*TY_qb12 - 2*TY_qc112)*TY_qc223)*
533                  pow(Z2,3))))*pow(Z1 + Z2,-1);
534       
535        TY_B21 = 6*phi*TY_qc221*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-(TY_qc221*(Z1 + Z2)*exp(Z1)*pow(Z2,2)) + 
536                  exp(2*Z2)*(exp(Z1)*(Z2*(2*TY_qb21*(-1 + Z2)*(Z1 + Z2) - Z2*(2*TY_qc121*Z2 + TY_qc221*(Z1 + Z2))) + TY_qa21*(Z1 + Z2)*(-2 + pow(Z2,2))) + 
537                  2*TY_qc121*(Z1 + Z2)*pow(Z2,2)) + 2*exp(Z2)*(-(TY_qc121*Z1*pow(Z2,2)) + (Z1 + Z2)*exp(Z1)*(TY_qa21 + (TY_qa21 + TY_qb21)*Z2 + TY_qc221*pow(Z2,2)))
538                  )*pow(Z1 + Z2,-1);
539       
540        TY_B22 = exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-12*phi*TY_qc221*TY_qc222*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 
541                  12*phi*exp(Z2)*((Z1 + Z2)*(TY_qa21*TY_qc222 + TY_qa22*TY_qc221*(1 + Z2) + Z2*(TY_qb22*TY_qc221 + TY_qc222*(TY_qa21 + TY_qb21 + 2*TY_qc221*Z2)))*exp(Z1) + 
542                  (-(TY_qc122*TY_qc221) - TY_qc121*TY_qc222)*Z1*pow(Z2,2)) + 
543                  exp(2*Z2)*(12*phi*(TY_qc122*TY_qc221 + TY_qc121*TY_qc222)*(Z1 + Z2)*pow(Z2,2) + 
544                  exp(Z1)*(6*phi*(2*Z2*(TY_qb22*TY_qc221*(-1 + Z2)*(Z1 + Z2) + TY_qb21*TY_qc222*(-1 + Z2)*(Z1 + Z2) - 
545                  Z2*(TY_qc221*TY_qc222*Z1 + TY_qc121*TY_qc222*Z2 + TY_qc221*(TY_qc122 + TY_qc222)*Z2)) + TY_qa22*TY_qc221*(Z1 + Z2)*(-2 + pow(Z2,2)) + 
546                  TY_qa21*TY_qc222*(Z1 + Z2)*(-2 + pow(Z2,2))) + TY_q22*TY_qc221*(Z1 + Z2)*pow(Z2,3))))*pow(Z1 + Z2,-1);
547       
548        TY_B23 = exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-6*phi*(Z1 + Z2)*exp(Z1)*(2*TY_qc221*TY_qc223 + 2*TY_qc212*TY_qc232 + pow(TY_qc222,2))*pow(Z2,2) + 
549                  12*phi*exp(Z2)*((Z1 + Z2)*exp(Z1)*(TY_qa23*TY_qc221 + TY_qa22*TY_qc222 + TY_qa21*TY_qc223 + TY_qa12*TY_qc232 + TY_qa32*TY_qc212*(1 + Z2) + 
550                  Z2*(TY_qb32*TY_qc212 + (TY_qa23 + TY_qb23)*TY_qc221 + (TY_qa22 + TY_qb22)*TY_qc222 + (TY_qa21 + TY_qb21)*TY_qc223 + (TY_qa12 + TY_qb12)*TY_qc232 + 
551                  Z2*(TY_q22*TY_qc221 + 2*TY_qc221*TY_qc223 + 2*TY_qc212*TY_qc232 + pow(TY_qc222,2)))) + 
552                  (-(TY_qc132*TY_qc212) - TY_qc123*TY_qc221 - TY_qc122*TY_qc222 - TY_qc121*TY_qc223 - TY_qc112*TY_qc232)*Z1*pow(Z2,2)) + 
553                  exp(2*Z2)*(12*phi*(TY_qc132*TY_qc212 + TY_qc123*TY_qc221 + TY_qc122*TY_qc222 + TY_qc121*(TY_q22 + TY_qc223) + TY_qc112*TY_qc232)*(Z1 + Z2)*pow(Z2,2) + 
554                  exp(Z1)*(TY_q22*TY_qc222*(Z1 + Z2)*pow(Z2,3) - 6*phi*
555                  (2*(TY_qa32*TY_qc212 + TY_qa23*TY_qc221 + TY_qa22*TY_qc222 + TY_qa21*(TY_q22 + TY_qc223) + TY_qa12*TY_qc232)*Z1 + 
556                  2*(TY_qa32*TY_qc212 + TY_qa23*TY_qc221 + TY_qa22*TY_qc222 + TY_qa21*(TY_q22 + TY_qc223) + TY_qa12*TY_qc232 + 
557                  (TY_qb32*TY_qc212 + TY_qb23*TY_qc221 + TY_qb22*TY_qc222 + TY_qb21*(TY_q22 + TY_qc223) + TY_qb12*TY_qc232)*Z1)*Z2 - 
558                  (-2*(TY_qb32*TY_qc212 + TY_qb23*TY_qc221 + TY_qb22*TY_qc222 + TY_qb21*(TY_q22 + TY_qc223) + TY_qb12*TY_qc232) + 
559                  (TY_q22*(TY_qa21 + 2*TY_qb21 - 2*TY_qc221) + (TY_qa22 + 2*TY_qb22 - TY_qc222)*TY_qc222 + TY_qc221*(TY_qa23 + 2*TY_qb23 - 2*TY_qc223) + TY_qa21*TY_qc223 + 
560                  2*TY_qb21*TY_qc223 + TY_qc212*(TY_qa32 + 2*TY_qb32 - 2*TY_qc232) + TY_qa12*TY_qc232 + 2*TY_qb12*TY_qc232)*Z1)*pow(Z2,2) - 
561                  (TY_q22*(TY_qa21 + 2*TY_qb21 - 2*TY_qc121 - 2*TY_qc212 - 2*TY_qc221) + (TY_qa22 + 2*TY_qb22 - 2*TY_qc122 - TY_qc222)*TY_qc222 + 
562                  TY_qc221*(TY_qa23 + 2*TY_qb23 - 2*TY_qc123 - 2*TY_qc223) + (TY_qa21 + 2*TY_qb21 - 2*TY_qc121)*TY_qc223 + 
563                  TY_qc212*(TY_qa32 + 2*TY_qb32 - 2*TY_qc132 - 2*TY_qc232) + (TY_qa12 + 2*TY_qb12 - 2*TY_qc112)*TY_qc232)*pow(Z2,3)))))*pow(Z1 + Z2,-1);
564       
565        TY_B24 = exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(12*phi*(Z1 + Z2)*
566                  (TY_qa22*TY_qc223 + TY_qa23*TY_qc222*(1 + Z2) + Z2*(TY_qb23*TY_qc222 + (TY_qa22 + TY_qb22)*TY_qc223 + TY_qc222*(TY_q22 + 2*TY_qc223)*Z2))*exp(Z1 + Z2) - 
567                  12*phi*TY_qc222*TY_qc223*(Z1 + Z2)*exp(Z1)*pow(Z2,2) - 12*phi*(TY_qc123*TY_qc222 + TY_qc122*TY_qc223)*Z1*exp(Z2)*pow(Z2,2) + 
568                  12*phi*(TY_qc123*TY_qc222 + TY_qc122*(TY_q22 + TY_qc223))*(Z1 + Z2)*exp(2*Z2)*pow(Z2,2) + 
569                  exp(Z1 + 2*Z2)*(6*phi*(2*Z2*(TY_qb23*TY_qc222*(-1 + Z2)*(Z1 + Z2) + TY_qb22*TY_qc223*(-1 + Z2)*(Z1 + Z2) - 
570                  Z2*(TY_qc222*TY_qc223*Z1 + TY_qc122*TY_qc223*Z2 + TY_qc222*(TY_qc123 + TY_qc223)*Z2)) + TY_qa23*TY_qc222*(Z1 + Z2)*(-2 + pow(Z2,2)) + 
571                  TY_qa22*TY_qc223*(Z1 + Z2)*(-2 + pow(Z2,2))) + (Z1 + Z2)*pow(TY_q22,2)*pow(Z2,3) + 
572                  TY_q22*(6*phi*(2*Z2*(TY_qb22*(-1 + Z2)*(Z1 + Z2) - Z2*(TY_qc222*Z1 + (TY_qc122 + TY_qc222)*Z2)) + TY_qa22*(Z1 + Z2)*(-2 + pow(Z2,2))) + 
573                  TY_qc223*(Z1 + Z2)*pow(Z2,3))))*pow(Z1 + Z2,-1);
574       
575        TY_B25 = -6*phi*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*((Z1 + Z2)*exp(Z1)*pow(TY_qc223,2)*pow(Z2,2) + 
576                  (TY_q22 + TY_qc223)*exp(2*Z2)*(exp(Z1)*(Z2*(-2*TY_qb23*(-1 + Z2)*(Z1 + Z2) + Z2*((TY_q22 + TY_qc223)*Z1 + (TY_q22 + 2*TY_qc123 + TY_qc223)*Z2)) - 
577                  TY_qa23*(Z1 + Z2)*(-2 + pow(Z2,2))) - 2*TY_qc123*(Z1 + Z2)*pow(Z2,2)) + 
578                  2*TY_qc223*exp(Z2)*(TY_qc123*Z1*pow(Z2,2) - (Z1 + Z2)*exp(Z1)*(TY_qa23 + (TY_qa23 + TY_qb23)*Z2 + (TY_q22 + TY_qc223)*pow(Z2,2))))*pow(Z1 + Z2,-1);
579       
580        TY_B32 = 6*phi*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-2*TY_qc221*TY_qc232*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 
581                  2*exp(Z2)*((Z1 + Z2)*(TY_qa21*TY_qc232 + TY_qa32*TY_qc221*(1 + Z2) + Z2*(TY_qb32*TY_qc221 + TY_qc232*(TY_qa21 + TY_qb21 + 2*TY_qc221*Z2)))*exp(Z1) + 
582                  (-(TY_qc132*TY_qc221) - TY_qc121*TY_qc232)*Z1*pow(Z2,2)) + 
583                  exp(2*Z2)*(exp(Z1)*(2*Z2*(TY_qb32*TY_qc221*(-1 + Z2)*(Z1 + Z2) + TY_qb21*TY_qc232*(-1 + Z2)*(Z1 + Z2) - 
584                  Z2*(TY_qc221*TY_qc232*Z1 + TY_qc121*TY_qc232*Z2 + TY_qc221*(TY_q22 + TY_qc132 + TY_qc232)*Z2)) + TY_qa32*TY_qc221*(Z1 + Z2)*(-2 + pow(Z2,2)) + 
585                  TY_qa21*TY_qc232*(Z1 + Z2)*(-2 + pow(Z2,2))) + 2*(TY_qc132*TY_qc221 + TY_qc121*TY_qc232)*(Z1 + Z2)*pow(Z2,2)))*pow(Z1 + Z2,-1);
586       
587        TY_B34 = -6*phi*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(2*TY_qc223*TY_qc232*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 
588                  2*exp(Z2)*(-((Z1 + Z2)*(TY_qa23*TY_qc232 + TY_qa32*TY_qc223*(1 + Z2) + Z2*(TY_qb32*TY_qc223 + TY_qc232*(TY_qa23 + TY_qb23 + TY_q22*Z2 + 2*TY_qc223*Z2)))*exp(Z1)) + 
589                  (TY_qc132*TY_qc223 + TY_qc123*TY_qc232)*Z1*pow(Z2,2)) + exp(2*Z2)*
590                  (-2*(TY_qc132*(TY_q22 + TY_qc223) + TY_qc123*TY_qc232)*(Z1 + Z2)*pow(Z2,2) + 
591                  exp(Z1)*(2*(TY_qa32*(TY_q22 + TY_qc223) + TY_qa23*TY_qc232)*Z1 + 
592                  2*(TY_qa32*(TY_q22 + TY_qc223) + TY_qa23*TY_qc232 + (TY_qb32*(TY_q22 + TY_qc223) + TY_qb23*TY_qc232)*Z1)*Z2 - 
593                  (-2*(TY_qb32*(TY_q22 + TY_qc223) + TY_qb23*TY_qc232) + ((TY_qa32 + 2*TY_qb32)*(TY_q22 + TY_qc223) + (-2*TY_q22 + TY_qa23 + 2*TY_qb23 - 2*TY_qc223)*TY_qc232)*Z1)*
594                  pow(Z2,2) + ((2*TY_q22 - TY_qa32 - 2*TY_qb32 + 2*TY_qc132)*(TY_q22 + TY_qc223) + (2*TY_q22 - TY_qa23 + 2*(-TY_qb23 + TY_qc123 + TY_qc223))*TY_qc232)*pow(Z2,3))))*pow(Z1 + Z2,-1);
595       
596        /*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));
597       
598        TY_B12 /= norm_B;
599        TY_B14 /= norm_B;
600        TY_B21 /= norm_B;
601        TY_B22 /= norm_B;
602        TY_B23 /= norm_B;
603        TY_B24 /= norm_B;
604        TY_B25 /= norm_B;
605        TY_B32 /= norm_B;
606        TY_B34 /= norm_B; */
607/*     
608        if( debug )
609        {
610                sprintf(buf,"\rNonlinear equation 2 = \r" );
611                XOPNotice(buf);
612                sprintf(buf, "%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0.,  TY_B12, 0.,  TY_B14, 0.  );
613                XOPNotice(buf);
614                sprintf(buf, "%f\t\t%f\t\t%f\t\t%f\t\t%f\r", TY_B21, TY_B22, TY_B23, TY_B24, TY_B25 );
615                XOPNotice(buf);
616                sprintf(buf, "%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0.,  TY_B32, 0.,  TY_B34, 0.  );
617                XOPNotice(buf);
618        }
619*/
620       
621        /* decrease order of nonlinear equation 1 by means of equation 2 */
622       
623        TY_F14 = -(TY_A32*TY_B12*TY_B32) + TY_A52*pow(TY_B12,2) + TY_A12*pow(TY_B32,2);
624        TY_F16 = 2*TY_A52*TY_B12*TY_B14 - TY_A32*TY_B14*TY_B32 - TY_A32*TY_B12*TY_B34 + 2*TY_A12*TY_B32*TY_B34;
625        TY_F18 = -(TY_A32*TY_B14*TY_B34) + TY_A52*pow(TY_B14,2) + TY_A12*pow(TY_B34,2);
626        TY_F23 = 2*TY_A52*TY_B12*TY_B21 - TY_A41*TY_B12*TY_B32 - TY_A32*TY_B21*TY_B32 + TY_A21*pow(TY_B32,2);
627        TY_F24 = 2*TY_A52*TY_B12*TY_B22 - TY_A42*TY_B12*TY_B32 - TY_A32*TY_B22*TY_B32 + TY_A22*pow(TY_B32,2);
628        TY_F25 = 2*TY_A52*TY_B14*TY_B21 + 2*TY_A52*TY_B12*TY_B23 - TY_A43*TY_B12*TY_B32 - TY_A41*TY_B14*TY_B32 - TY_A32*TY_B23*TY_B32 - TY_A41*TY_B12*TY_B34 - TY_A32*TY_B21*TY_B34 + 2*TY_A21*TY_B32*TY_B34 + TY_A23*pow(TY_B32,2);
629        TY_F26 = 2*TY_A52*TY_B14*TY_B22 + 2*TY_A52*TY_B12*TY_B24 - TY_A42*TY_B14*TY_B32 - TY_A32*TY_B24*TY_B32 - TY_A42*TY_B12*TY_B34 - TY_A32*TY_B22*TY_B34 + 2*TY_A22*TY_B32*TY_B34;
630        TY_F27 = 2*TY_A52*TY_B14*TY_B23 + 2*TY_A52*TY_B12*TY_B25 - TY_A43*TY_B14*TY_B32 - TY_A32*TY_B25*TY_B32 - TY_A43*TY_B12*TY_B34 - TY_A41*TY_B14*TY_B34 - TY_A32*TY_B23*TY_B34 + 2*TY_A23*TY_B32*TY_B34 + TY_A21*pow(TY_B34,2);
631        TY_F28 = 2*TY_A52*TY_B14*TY_B24 - TY_A42*TY_B14*TY_B34 - TY_A32*TY_B24*TY_B34 + TY_A22*pow(TY_B34,2);
632        TY_F29 = 2*TY_A52*TY_B14*TY_B25 - TY_A43*TY_B14*TY_B34 - TY_A32*TY_B25*TY_B34 + TY_A23*pow(TY_B34,2);
633        TY_F32 = -(TY_A41*TY_B21*TY_B32) + TY_A52*pow(TY_B21,2);
634        TY_F33 = 2*TY_A52*TY_B21*TY_B22 - TY_A42*TY_B21*TY_B32 - TY_A41*TY_B22*TY_B32;
635        TY_F34 = 2*TY_A52*TY_B21*TY_B23 - TY_A43*TY_B21*TY_B32 - TY_A42*TY_B22*TY_B32 - TY_A41*TY_B23*TY_B32 - TY_A41*TY_B21*TY_B34 + TY_A52*pow(TY_B22,2);
636        TY_F35 = 2*TY_A52*TY_B22*TY_B23 + 2*TY_A52*TY_B21*TY_B24 - TY_A43*TY_B22*TY_B32 - TY_A42*TY_B23*TY_B32 - TY_A41*TY_B24*TY_B32 - TY_A42*TY_B21*TY_B34 - TY_A41*TY_B22*TY_B34;
637        TY_F36 = 2*TY_A52*TY_B22*TY_B24 + 2*TY_A52*TY_B21*TY_B25 - TY_A43*TY_B23*TY_B32 - TY_A42*TY_B24*TY_B32 - TY_A41*TY_B25*TY_B32 - TY_A43*TY_B21*TY_B34 - TY_A42*TY_B22*TY_B34 - TY_A41*TY_B23*TY_B34 + TY_A52*pow(TY_B23,2);
638        TY_F37 = 2*TY_A52*TY_B23*TY_B24 + 2*TY_A52*TY_B22*TY_B25 - TY_A43*TY_B24*TY_B32 - TY_A42*TY_B25*TY_B32 - TY_A43*TY_B22*TY_B34 - TY_A42*TY_B23*TY_B34 - TY_A41*TY_B24*TY_B34;
639        TY_F38 = 2*TY_A52*TY_B23*TY_B25 - TY_A43*TY_B25*TY_B32 - TY_A43*TY_B23*TY_B34 - TY_A42*TY_B24*TY_B34 - TY_A41*TY_B25*TY_B34 + TY_A52*pow(TY_B24,2);
640        TY_F39 = 2*TY_A52*TY_B24*TY_B25 - TY_A43*TY_B24*TY_B34 - TY_A42*TY_B25*TY_B34;
641        TY_F310 = -(TY_A43*TY_B25*TY_B34) + TY_A52*pow(TY_B25,2);
642
643/*
644        if( debug )
645        {
646                sprintf(buf,"\rF = \r" );
647                XOPNotice(buf);
648                sprintf(buf, "%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.    );
649                XOPNotice(buf);
650                sprintf(buf, "%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.    );
651                XOPNotice(buf);
652                sprintf(buf, "%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 );
653                XOPNotice(buf);
654        }
655*/
656 
657        TY_G13  = -(TY_B12*TY_F32);
658        TY_G14  = -(TY_B12*TY_F33);
659        TY_G15  = TY_B32*TY_F14 - TY_B14*TY_F32 - TY_B12*TY_F34;
660        TY_G16  = -(TY_B14*TY_F33) - TY_B12*TY_F35;
661        TY_G17  = TY_B34*TY_F14 + TY_B32*TY_F16 - TY_B14*TY_F34 - TY_B12*TY_F36;
662        TY_G18  = -(TY_B14*TY_F35) - TY_B12*TY_F37;
663        TY_G19  = TY_B34*TY_F16 + TY_B32*TY_F18 - TY_B14*TY_F36 - TY_B12*TY_F38;
664        TY_G110 = -(TY_B14*TY_F37) - TY_B12*TY_F39;
665        TY_G111 = TY_B34*TY_F18 - TY_B12*TY_F310 - TY_B14*TY_F38;
666        TY_G112 = -(TY_B14*TY_F39);
667        TY_G113 = -(TY_B14*TY_F310);
668        TY_G22  = -(TY_B21*TY_F32);
669        TY_G23  = -(TY_B22*TY_F32) - TY_B21*TY_F33;
670        TY_G24  = TY_B32*TY_F23 - TY_B23*TY_F32 - TY_B22*TY_F33 - TY_B21*TY_F34;
671        TY_G25  = TY_B32*TY_F24 - TY_B24*TY_F32 - TY_B23*TY_F33 - TY_B22*TY_F34 - TY_B21*TY_F35;
672        TY_G26  = TY_B34*TY_F23 + TY_B32*TY_F25 - TY_B25*TY_F32 - TY_B24*TY_F33 - TY_B23*TY_F34 - TY_B22*TY_F35 - TY_B21*TY_F36;
673        TY_G27  = TY_B34*TY_F24 + TY_B32*TY_F26 - TY_B25*TY_F33 - TY_B24*TY_F34 - TY_B23*TY_F35 - TY_B22*TY_F36 - TY_B21*TY_F37;
674        TY_G28  = TY_B34*TY_F25 + TY_B32*TY_F27 - TY_B25*TY_F34 - TY_B24*TY_F35 - TY_B23*TY_F36 - TY_B22*TY_F37 - TY_B21*TY_F38;
675        TY_G29  = TY_B34*TY_F26 + TY_B32*TY_F28 - TY_B25*TY_F35 - TY_B24*TY_F36 - TY_B23*TY_F37 - TY_B22*TY_F38 - TY_B21*TY_F39;
676        TY_G210 = TY_B34*TY_F27 + TY_B32*TY_F29 - TY_B21*TY_F310 - TY_B25*TY_F36 - TY_B24*TY_F37 - TY_B23*TY_F38 - TY_B22*TY_F39;
677        TY_G211 = TY_B34*TY_F28 - TY_B22*TY_F310 - TY_B25*TY_F37 - TY_B24*TY_F38 - TY_B23*TY_F39;
678        TY_G212 = TY_B34*TY_F29 - TY_B23*TY_F310 - TY_B25*TY_F38 - TY_B24*TY_F39;
679        TY_G213 = -(TY_B24*TY_F310) - TY_B25*TY_F39;
680        TY_G214 = -(TY_B25*TY_F310);
681/*     
682        if( debug )
683        {
684                sprintf(buf,"\rG = \r" );
685                XOPNotice(buf);
686                sprintf(buf, "%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. );
687                XOPNotice(buf);
688                sprintf(buf, "%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 );
689                XOPNotice(buf);
690        }
691*/     
692        // coefficients for polynomial
693        TY_w[0] = (-(TY_A21*TY_B12) + TY_A12*TY_B21)*(TY_A52*TY_B21 - TY_A41*TY_B32)*pow(TY_B21,2)*pow(TY_B32,3);
694       
695        TY_w[1] = 2*TY_B32*TY_G13*TY_G14 - TY_B24*TY_G13*TY_G22 - TY_B23*TY_G14*TY_G22 - TY_B22*TY_G15*TY_G22 - TY_B21*TY_G16*TY_G22 - TY_B23*TY_G13*TY_G23 - TY_B22*TY_G14*TY_G23 - TY_B21*TY_G15*TY_G23 + 2*TY_B14*TY_G22*TY_G23 - 
696        TY_B22*TY_G13*TY_G24 - TY_B21*TY_G14*TY_G24 + 2*TY_B12*TY_G23*TY_G24 - TY_B21*TY_G13*TY_G25 + 2*TY_B12*TY_G22*TY_G25;
697       
698        TY_w[2] = -(TY_B25*TY_G13*TY_G22) - TY_B24*TY_G14*TY_G22 - TY_B23*TY_G15*TY_G22 - TY_B22*TY_G16*TY_G22 - TY_B21*TY_G17*TY_G22 - TY_B24*TY_G13*TY_G23 - TY_B23*TY_G14*TY_G23 - TY_B22*TY_G15*TY_G23 - TY_B21*TY_G16*TY_G23 - 
699        TY_B23*TY_G13*TY_G24 - TY_B22*TY_G14*TY_G24 - TY_B21*TY_G15*TY_G24 + 2*TY_B14*TY_G22*TY_G24 - TY_B22*TY_G13*TY_G25 - TY_B21*TY_G14*TY_G25 + 2*TY_B12*TY_G23*TY_G25 - TY_B21*TY_G13*TY_G26 + 2*TY_B12*TY_G22*TY_G26 + 
700        TY_B34*pow(TY_G13,2) + TY_B32*(2*TY_G13*TY_G15 + pow(TY_G14,2)) + TY_B14*pow(TY_G23,2) + TY_B12*pow(TY_G24,2);
701       
702        TY_w[3] = 2*TY_B34*TY_G13*TY_G14 + 2*TY_B32*(TY_G14*TY_G15 + TY_G13*TY_G16) - TY_B25*TY_G14*TY_G22 - TY_B24*TY_G15*TY_G22 - TY_B23*TY_G16*TY_G22 - TY_B22*TY_G17*TY_G22 - TY_B21*TY_G18*TY_G22 - TY_B25*TY_G13*TY_G23 - 
703        TY_B24*TY_G14*TY_G23 - TY_B23*TY_G15*TY_G23 - TY_B22*TY_G16*TY_G23 - TY_B21*TY_G17*TY_G23 - TY_B24*TY_G13*TY_G24 - TY_B23*TY_G14*TY_G24 - TY_B22*TY_G15*TY_G24 - TY_B21*TY_G16*TY_G24 + 2*TY_B14*TY_G23*TY_G24 - 
704        TY_B23*TY_G13*TY_G25 - TY_B22*TY_G14*TY_G25 - TY_B21*TY_G15*TY_G25 + 2*TY_B14*TY_G22*TY_G25 + 2*TY_B12*TY_G24*TY_G25 - TY_B22*TY_G13*TY_G26 - TY_B21*TY_G14*TY_G26 + 2*TY_B12*TY_G23*TY_G26 - TY_B21*TY_G13*TY_G27 + 
705        2*TY_B12*TY_G22*TY_G27;
706       
707        TY_w[4] = -(TY_B25*TY_G15*TY_G22) - TY_B24*TY_G16*TY_G22 - TY_B23*TY_G17*TY_G22 - TY_B22*TY_G18*TY_G22 - TY_B21*TY_G19*TY_G22 - TY_B25*TY_G14*TY_G23 - TY_B24*TY_G15*TY_G23 - TY_B23*TY_G16*TY_G23 - TY_B22*TY_G17*TY_G23 - 
708        TY_B21*TY_G18*TY_G23 - TY_B25*TY_G13*TY_G24 - TY_B24*TY_G14*TY_G24 - TY_B23*TY_G15*TY_G24 - TY_B22*TY_G16*TY_G24 - TY_B21*TY_G17*TY_G24 - TY_B24*TY_G13*TY_G25 - TY_B23*TY_G14*TY_G25 - TY_B22*TY_G15*TY_G25 - 
709        TY_B21*TY_G16*TY_G25 + 2*TY_B14*TY_G23*TY_G25 - TY_B23*TY_G13*TY_G26 - TY_B22*TY_G14*TY_G26 - TY_B21*TY_G15*TY_G26 + 2*TY_B14*TY_G22*TY_G26 + 2*TY_B12*TY_G24*TY_G26 - TY_B22*TY_G13*TY_G27 - TY_B21*TY_G14*TY_G27 + 
710        2*TY_B12*TY_G23*TY_G27 - TY_B21*TY_G13*TY_G28 + 2*TY_B12*TY_G22*TY_G28 + TY_B34*(2*TY_G13*TY_G15 + pow(TY_G14,2)) + TY_B32*(2*TY_G14*TY_G16 + 2*TY_G13*TY_G17 + pow(TY_G15,2)) + TY_B14*pow(TY_G24,2) + 
711        TY_B12*pow(TY_G25,2);
712       
713        TY_w[5] = 2*TY_B34*(TY_G14*TY_G15 + TY_G13*TY_G16) + 2*TY_B32*(TY_G15*TY_G16 + TY_G14*TY_G17 + TY_G13*TY_G18) - TY_B21*TY_G110*TY_G22 - TY_B25*TY_G16*TY_G22 - TY_B24*TY_G17*TY_G22 - TY_B23*TY_G18*TY_G22 - TY_B22*TY_G19*TY_G22 - 
714        TY_B25*TY_G15*TY_G23 - TY_B24*TY_G16*TY_G23 - TY_B23*TY_G17*TY_G23 - TY_B22*TY_G18*TY_G23 - TY_B21*TY_G19*TY_G23 - TY_B25*TY_G14*TY_G24 - TY_B24*TY_G15*TY_G24 - TY_B23*TY_G16*TY_G24 - TY_B22*TY_G17*TY_G24 - 
715        TY_B21*TY_G18*TY_G24 - TY_B25*TY_G13*TY_G25 - TY_B24*TY_G14*TY_G25 - TY_B23*TY_G15*TY_G25 - TY_B22*TY_G16*TY_G25 - TY_B21*TY_G17*TY_G25 + 2*TY_B14*TY_G24*TY_G25 - TY_B24*TY_G13*TY_G26 - TY_B23*TY_G14*TY_G26 - 
716        TY_B22*TY_G15*TY_G26 - TY_B21*TY_G16*TY_G26 + 2*TY_B14*TY_G23*TY_G26 + 2*TY_B12*TY_G25*TY_G26 - TY_B23*TY_G13*TY_G27 - TY_B22*TY_G14*TY_G27 - TY_B21*TY_G15*TY_G27 + 2*TY_B14*TY_G22*TY_G27 + 2*TY_B12*TY_G24*TY_G27 - 
717        TY_B22*TY_G13*TY_G28 - TY_B21*TY_G14*TY_G28 + 2*TY_B12*TY_G23*TY_G28 - TY_B21*TY_G13*TY_G29 + 2*TY_B12*TY_G22*TY_G29;
718       
719        TY_w[6] = -(TY_B22*TY_G110*TY_G22) - TY_B21*TY_G111*TY_G22 - TY_B25*TY_G17*TY_G22 - TY_B24*TY_G18*TY_G22 - TY_B23*TY_G19*TY_G22 + TY_G210*(-(TY_B21*TY_G13) + 2*TY_B12*TY_G22) - TY_B21*TY_G110*TY_G23 - TY_B25*TY_G16*TY_G23 - 
720        TY_B24*TY_G17*TY_G23 - TY_B23*TY_G18*TY_G23 - TY_B22*TY_G19*TY_G23 - TY_B25*TY_G15*TY_G24 - TY_B24*TY_G16*TY_G24 - TY_B23*TY_G17*TY_G24 - TY_B22*TY_G18*TY_G24 - TY_B21*TY_G19*TY_G24 - TY_B25*TY_G14*TY_G25 - 
721        TY_B24*TY_G15*TY_G25 - TY_B23*TY_G16*TY_G25 - TY_B22*TY_G17*TY_G25 - TY_B21*TY_G18*TY_G25 - TY_B25*TY_G13*TY_G26 - TY_B24*TY_G14*TY_G26 - TY_B23*TY_G15*TY_G26 - TY_B22*TY_G16*TY_G26 - TY_B21*TY_G17*TY_G26 + 
722        2*TY_B14*TY_G24*TY_G26 - TY_B24*TY_G13*TY_G27 - TY_B23*TY_G14*TY_G27 - TY_B22*TY_G15*TY_G27 - TY_B21*TY_G16*TY_G27 + 2*TY_B14*TY_G23*TY_G27 + 2*TY_B12*TY_G25*TY_G27 - TY_B23*TY_G13*TY_G28 - TY_B22*TY_G14*TY_G28 - 
723        TY_B21*TY_G15*TY_G28 + 2*TY_B14*TY_G22*TY_G28 + 2*TY_B12*TY_G24*TY_G28 - TY_B22*TY_G13*TY_G29 - TY_B21*TY_G14*TY_G29 + 2*TY_B12*TY_G23*TY_G29 + TY_B34*(2*TY_G14*TY_G16 + 2*TY_G13*TY_G17 + pow(TY_G15,2)) + 
724        TY_B32*(2*(TY_G15*TY_G17 + TY_G14*TY_G18 + TY_G13*TY_G19) + pow(TY_G16,2)) + TY_B14*pow(TY_G25,2) + TY_B12*pow(TY_G26,2);
725       
726        TY_w[7] = 2*TY_B34*(TY_G15*TY_G16 + TY_G14*TY_G17 + TY_G13*TY_G18) + 2*TY_B32*(TY_G110*TY_G13 + TY_G16*TY_G17 + TY_G15*TY_G18 + TY_G14*TY_G19) - TY_B22*TY_G13*TY_G210 - TY_B21*TY_G14*TY_G210 - TY_B23*TY_G110*TY_G22 - 
727        TY_B22*TY_G111*TY_G22 - TY_B21*TY_G112*TY_G22 - TY_B25*TY_G18*TY_G22 - TY_B24*TY_G19*TY_G22 + TY_G211*(-(TY_B21*TY_G13) + 2*TY_B12*TY_G22) - TY_B22*TY_G110*TY_G23 - TY_B21*TY_G111*TY_G23 - TY_B25*TY_G17*TY_G23 - 
728        TY_B24*TY_G18*TY_G23 - TY_B23*TY_G19*TY_G23 + 2*TY_B12*TY_G210*TY_G23 - TY_B21*TY_G110*TY_G24 - TY_B25*TY_G16*TY_G24 - TY_B24*TY_G17*TY_G24 - TY_B23*TY_G18*TY_G24 - TY_B22*TY_G19*TY_G24 - TY_B25*TY_G15*TY_G25 - 
729        TY_B24*TY_G16*TY_G25 - TY_B23*TY_G17*TY_G25 - TY_B22*TY_G18*TY_G25 - TY_B21*TY_G19*TY_G25 - TY_B25*TY_G14*TY_G26 - TY_B24*TY_G15*TY_G26 - TY_B23*TY_G16*TY_G26 - TY_B22*TY_G17*TY_G26 - TY_B21*TY_G18*TY_G26 + 
730        2*TY_B14*TY_G25*TY_G26 - TY_B25*TY_G13*TY_G27 - TY_B24*TY_G14*TY_G27 - TY_B23*TY_G15*TY_G27 - TY_B22*TY_G16*TY_G27 - TY_B21*TY_G17*TY_G27 + 2*TY_B14*TY_G24*TY_G27 + 2*TY_B12*TY_G26*TY_G27 - TY_B24*TY_G13*TY_G28 - 
731        TY_B23*TY_G14*TY_G28 - TY_B22*TY_G15*TY_G28 - TY_B21*TY_G16*TY_G28 + 2*TY_B14*TY_G23*TY_G28 + 2*TY_B12*TY_G25*TY_G28 - TY_B23*TY_G13*TY_G29 - TY_B22*TY_G14*TY_G29 - TY_B21*TY_G15*TY_G29 + 2*TY_B14*TY_G22*TY_G29 + 
732        2*TY_B12*TY_G24*TY_G29;
733       
734        TY_w[8] = -(TY_B23*TY_G13*TY_G210) - TY_B22*TY_G14*TY_G210 - TY_B21*TY_G15*TY_G210 - TY_B22*TY_G13*TY_G211 - TY_B21*TY_G14*TY_G211 - TY_B21*TY_G13*TY_G212 - TY_B24*TY_G110*TY_G22 - TY_B23*TY_G111*TY_G22 - TY_B22*TY_G112*TY_G22 - 
735        TY_B21*TY_G113*TY_G22 - TY_B25*TY_G19*TY_G22 + 2*TY_B14*TY_G210*TY_G22 + 2*TY_B12*TY_G212*TY_G22 - TY_B23*TY_G110*TY_G23 - TY_B22*TY_G111*TY_G23 - TY_B21*TY_G112*TY_G23 - TY_B25*TY_G18*TY_G23 - TY_B24*TY_G19*TY_G23 + 
736        2*TY_B12*TY_G211*TY_G23 - TY_B22*TY_G110*TY_G24 - TY_B21*TY_G111*TY_G24 - TY_B25*TY_G17*TY_G24 - TY_B24*TY_G18*TY_G24 - TY_B23*TY_G19*TY_G24 + 2*TY_B12*TY_G210*TY_G24 - TY_B21*TY_G110*TY_G25 - TY_B25*TY_G16*TY_G25 - 
737        TY_B24*TY_G17*TY_G25 - TY_B23*TY_G18*TY_G25 - TY_B22*TY_G19*TY_G25 - TY_B25*TY_G15*TY_G26 - TY_B24*TY_G16*TY_G26 - TY_B23*TY_G17*TY_G26 - TY_B22*TY_G18*TY_G26 - TY_B21*TY_G19*TY_G26 - TY_B25*TY_G14*TY_G27 - 
738        TY_B24*TY_G15*TY_G27 - TY_B23*TY_G16*TY_G27 - TY_B22*TY_G17*TY_G27 - TY_B21*TY_G18*TY_G27 + 2*TY_B14*TY_G25*TY_G27 - TY_B25*TY_G13*TY_G28 - TY_B24*TY_G14*TY_G28 - TY_B23*TY_G15*TY_G28 - TY_B22*TY_G16*TY_G28 - 
739        TY_B21*TY_G17*TY_G28 + 2*TY_B14*TY_G24*TY_G28 + 2*TY_B12*TY_G26*TY_G28 - TY_B24*TY_G13*TY_G29 - TY_B23*TY_G14*TY_G29 - TY_B22*TY_G15*TY_G29 - TY_B21*TY_G16*TY_G29 + 2*TY_B14*TY_G23*TY_G29 + 2*TY_B12*TY_G25*TY_G29 + 
740        TY_B34*(2*(TY_G15*TY_G17 + TY_G14*TY_G18 + TY_G13*TY_G19) + pow(TY_G16,2)) + TY_B32*(2*(TY_G111*TY_G13 + TY_G110*TY_G14 + TY_G16*TY_G18 + TY_G15*TY_G19) + pow(TY_G17,2)) + TY_B14*pow(TY_G26,2) + 
741        TY_B12*pow(TY_G27,2);
742       
743        TY_w[9] = 2*TY_B34*(TY_G110*TY_G13 + TY_G16*TY_G17 + TY_G15*TY_G18 + TY_G14*TY_G19) + 2*TY_B32*(TY_G112*TY_G13 + TY_G111*TY_G14 + TY_G110*TY_G15 + TY_G17*TY_G18 + TY_G16*TY_G19) - TY_B24*TY_G13*TY_G210 - TY_B23*TY_G14*TY_G210 - 
744        TY_B22*TY_G15*TY_G210 - TY_B21*TY_G16*TY_G210 - TY_B23*TY_G13*TY_G211 - TY_B22*TY_G14*TY_G211 - TY_B21*TY_G15*TY_G211 - TY_B22*TY_G13*TY_G212 - TY_B21*TY_G14*TY_G212 - TY_B25*TY_G110*TY_G22 - TY_B24*TY_G111*TY_G22 - 
745        TY_B23*TY_G112*TY_G22 - TY_B22*TY_G113*TY_G22 + 2*TY_B14*TY_G211*TY_G22 + TY_G213*(-(TY_B21*TY_G13) + 2*TY_B12*TY_G22) - TY_B24*TY_G110*TY_G23 - TY_B23*TY_G111*TY_G23 - TY_B22*TY_G112*TY_G23 - 
746        TY_B21*TY_G113*TY_G23 - TY_B25*TY_G19*TY_G23 + 2*TY_B14*TY_G210*TY_G23 + 2*TY_B12*TY_G212*TY_G23 - TY_B23*TY_G110*TY_G24 - TY_B22*TY_G111*TY_G24 - TY_B21*TY_G112*TY_G24 - TY_B25*TY_G18*TY_G24 - TY_B24*TY_G19*TY_G24 + 
747        2*TY_B12*TY_G211*TY_G24 - TY_B22*TY_G110*TY_G25 - TY_B21*TY_G111*TY_G25 - TY_B25*TY_G17*TY_G25 - TY_B24*TY_G18*TY_G25 - TY_B23*TY_G19*TY_G25 + 2*TY_B12*TY_G210*TY_G25 - TY_B21*TY_G110*TY_G26 - TY_B25*TY_G16*TY_G26 - 
748        TY_B24*TY_G17*TY_G26 - TY_B23*TY_G18*TY_G26 - TY_B22*TY_G19*TY_G26 - TY_B25*TY_G15*TY_G27 - TY_B24*TY_G16*TY_G27 - TY_B23*TY_G17*TY_G27 - TY_B22*TY_G18*TY_G27 - TY_B21*TY_G19*TY_G27 + 2*TY_B14*TY_G26*TY_G27 - 
749        TY_B25*TY_G14*TY_G28 - TY_B24*TY_G15*TY_G28 - TY_B23*TY_G16*TY_G28 - TY_B22*TY_G17*TY_G28 - TY_B21*TY_G18*TY_G28 + 2*TY_B14*TY_G25*TY_G28 + 2*TY_B12*TY_G27*TY_G28 - TY_B25*TY_G13*TY_G29 - TY_B24*TY_G14*TY_G29 - 
750        TY_B23*TY_G15*TY_G29 - TY_B22*TY_G16*TY_G29 - TY_B21*TY_G17*TY_G29 + 2*TY_B14*TY_G24*TY_G29 + 2*TY_B12*TY_G26*TY_G29;
751       
752        TY_w[10] = -(TY_B25*TY_G13*TY_G210) - TY_B24*TY_G14*TY_G210 - TY_B23*TY_G15*TY_G210 - TY_B22*TY_G16*TY_G210 - TY_B21*TY_G17*TY_G210 - TY_B24*TY_G13*TY_G211 - TY_B23*TY_G14*TY_G211 - TY_B22*TY_G15*TY_G211 - TY_B21*TY_G16*TY_G211 - 
753        TY_B23*TY_G13*TY_G212 - TY_B22*TY_G14*TY_G212 - TY_B21*TY_G15*TY_G212 - TY_B22*TY_G13*TY_G213 - TY_B21*TY_G14*TY_G213 - TY_B21*TY_G13*TY_G214 - TY_B25*TY_G111*TY_G22 - TY_B24*TY_G112*TY_G22 - TY_B23*TY_G113*TY_G22 + 
754        2*TY_B14*TY_G212*TY_G22 + 2*TY_B12*TY_G214*TY_G22 - TY_B25*TY_G110*TY_G23 - TY_B24*TY_G111*TY_G23 - TY_B23*TY_G112*TY_G23 - TY_B22*TY_G113*TY_G23 + 2*TY_B14*TY_G211*TY_G23 + 2*TY_B12*TY_G213*TY_G23 - 
755        TY_B24*TY_G110*TY_G24 - TY_B23*TY_G111*TY_G24 - TY_B22*TY_G112*TY_G24 - TY_B21*TY_G113*TY_G24 - TY_B25*TY_G19*TY_G24 + 2*TY_B14*TY_G210*TY_G24 + 2*TY_B12*TY_G212*TY_G24 - TY_B23*TY_G110*TY_G25 - 
756        TY_B22*TY_G111*TY_G25 - TY_B21*TY_G112*TY_G25 - TY_B25*TY_G18*TY_G25 - TY_B24*TY_G19*TY_G25 + 2*TY_B12*TY_G211*TY_G25 - TY_B22*TY_G110*TY_G26 - TY_B21*TY_G111*TY_G26 - TY_B25*TY_G17*TY_G26 - TY_B24*TY_G18*TY_G26 - 
757        TY_B23*TY_G19*TY_G26 + 2*TY_B12*TY_G210*TY_G26 - TY_B21*TY_G110*TY_G27 - TY_B25*TY_G16*TY_G27 - TY_B24*TY_G17*TY_G27 - TY_B23*TY_G18*TY_G27 - TY_B22*TY_G19*TY_G27 - TY_B25*TY_G15*TY_G28 - TY_B24*TY_G16*TY_G28 - 
758        TY_B23*TY_G17*TY_G28 - TY_B22*TY_G18*TY_G28 - TY_B21*TY_G19*TY_G28 + 2*TY_B14*TY_G26*TY_G28 - TY_B25*TY_G14*TY_G29 - TY_B24*TY_G15*TY_G29 - TY_B23*TY_G16*TY_G29 - TY_B22*TY_G17*TY_G29 - TY_B21*TY_G18*TY_G29 + 
759        2*TY_B14*TY_G25*TY_G29 + 2*TY_B12*TY_G27*TY_G29 + TY_B34*(2*(TY_G111*TY_G13 + TY_G110*TY_G14 + TY_G16*TY_G18 + TY_G15*TY_G19) + pow(TY_G17,2)) + 
760        TY_B32*(2*(TY_G113*TY_G13 + TY_G112*TY_G14 + TY_G111*TY_G15 + TY_G110*TY_G16 + TY_G17*TY_G19) + pow(TY_G18,2)) + TY_B14*pow(TY_G27,2) + TY_B12*pow(TY_G28,2);
761       
762        TY_w[11] = 2*TY_B34*(TY_G112*TY_G13 + TY_G111*TY_G14 + TY_G110*TY_G15 + TY_G17*TY_G18 + TY_G16*TY_G19) + 2*TY_B32*(TY_G113*TY_G14 + TY_G112*TY_G15 + TY_G111*TY_G16 + TY_G110*TY_G17 + TY_G18*TY_G19) - TY_B25*TY_G14*TY_G210 - 
763        TY_B24*TY_G15*TY_G210 - TY_B23*TY_G16*TY_G210 - TY_B22*TY_G17*TY_G210 - TY_B21*TY_G18*TY_G210 - TY_B25*TY_G13*TY_G211 - TY_B24*TY_G14*TY_G211 - TY_B23*TY_G15*TY_G211 - TY_B22*TY_G16*TY_G211 - TY_B21*TY_G17*TY_G211 - 
764        TY_B24*TY_G13*TY_G212 - TY_B23*TY_G14*TY_G212 - TY_B22*TY_G15*TY_G212 - TY_B21*TY_G16*TY_G212 - TY_B23*TY_G13*TY_G213 - TY_B22*TY_G14*TY_G213 - TY_B21*TY_G15*TY_G213 - TY_B25*TY_G112*TY_G22 - TY_B24*TY_G113*TY_G22 + 
765        2*TY_B14*TY_G213*TY_G22 - TY_B25*TY_G111*TY_G23 - TY_B24*TY_G112*TY_G23 - TY_B23*TY_G113*TY_G23 + 2*TY_B14*TY_G212*TY_G23 - TY_G214*(TY_B22*TY_G13 + TY_B21*TY_G14 - 2*TY_B12*TY_G23) - TY_B25*TY_G110*TY_G24 - 
766        TY_B24*TY_G111*TY_G24 - TY_B23*TY_G112*TY_G24 - TY_B22*TY_G113*TY_G24 + 2*TY_B14*TY_G211*TY_G24 + 2*TY_B12*TY_G213*TY_G24 - TY_B24*TY_G110*TY_G25 - TY_B23*TY_G111*TY_G25 - TY_B22*TY_G112*TY_G25 - 
767        TY_B21*TY_G113*TY_G25 - TY_B25*TY_G19*TY_G25 + 2*TY_B14*TY_G210*TY_G25 + 2*TY_B12*TY_G212*TY_G25 - TY_B23*TY_G110*TY_G26 - TY_B22*TY_G111*TY_G26 - TY_B21*TY_G112*TY_G26 - TY_B25*TY_G18*TY_G26 - TY_B24*TY_G19*TY_G26 + 
768        2*TY_B12*TY_G211*TY_G26 - TY_B22*TY_G110*TY_G27 - TY_B21*TY_G111*TY_G27 - TY_B25*TY_G17*TY_G27 - TY_B24*TY_G18*TY_G27 - TY_B23*TY_G19*TY_G27 + 2*TY_B12*TY_G210*TY_G27 - TY_B21*TY_G110*TY_G28 - TY_B25*TY_G16*TY_G28 - 
769        TY_B24*TY_G17*TY_G28 - TY_B23*TY_G18*TY_G28 - TY_B22*TY_G19*TY_G28 + 2*TY_B14*TY_G27*TY_G28 - TY_B25*TY_G15*TY_G29 - TY_B24*TY_G16*TY_G29 - TY_B23*TY_G17*TY_G29 - TY_B22*TY_G18*TY_G29 - TY_B21*TY_G19*TY_G29 + 
770        2*TY_B14*TY_G26*TY_G29 + 2*TY_B12*TY_G28*TY_G29;
771       
772        TY_w[12] = -(TY_B25*TY_G15*TY_G210) - TY_B24*TY_G16*TY_G210 - TY_B23*TY_G17*TY_G210 - TY_B22*TY_G18*TY_G210 - TY_B21*TY_G19*TY_G210 - TY_B25*TY_G14*TY_G211 - TY_B24*TY_G15*TY_G211 - TY_B23*TY_G16*TY_G211 - TY_B22*TY_G17*TY_G211 - 
773        TY_B21*TY_G18*TY_G211 - TY_B25*TY_G13*TY_G212 - TY_B24*TY_G14*TY_G212 - TY_B23*TY_G15*TY_G212 - TY_B22*TY_G16*TY_G212 - TY_B21*TY_G17*TY_G212 - TY_B24*TY_G13*TY_G213 - TY_B23*TY_G14*TY_G213 - TY_B22*TY_G15*TY_G213 - 
774        TY_B21*TY_G16*TY_G213 - TY_B25*TY_G113*TY_G22 - TY_B25*TY_G112*TY_G23 - TY_B24*TY_G113*TY_G23 + 2*TY_B14*TY_G213*TY_G23 - TY_B25*TY_G111*TY_G24 - TY_B24*TY_G112*TY_G24 - TY_B23*TY_G113*TY_G24 + 
775        2*TY_B14*TY_G212*TY_G24 - TY_G214*(TY_B23*TY_G13 + TY_B22*TY_G14 + TY_B21*TY_G15 - 2*TY_B14*TY_G22 - 2*TY_B12*TY_G24) - TY_B25*TY_G110*TY_G25 - TY_B24*TY_G111*TY_G25 - TY_B23*TY_G112*TY_G25 - 
776        TY_B22*TY_G113*TY_G25 + 2*TY_B14*TY_G211*TY_G25 + 2*TY_B12*TY_G213*TY_G25 - TY_B24*TY_G110*TY_G26 - TY_B23*TY_G111*TY_G26 - TY_B22*TY_G112*TY_G26 - TY_B21*TY_G113*TY_G26 - TY_B25*TY_G19*TY_G26 + 
777        2*TY_B14*TY_G210*TY_G26 + 2*TY_B12*TY_G212*TY_G26 - TY_B23*TY_G110*TY_G27 - TY_B22*TY_G111*TY_G27 - TY_B21*TY_G112*TY_G27 - TY_B25*TY_G18*TY_G27 - TY_B24*TY_G19*TY_G27 + 2*TY_B12*TY_G211*TY_G27 - 
778        TY_B22*TY_G110*TY_G28 - TY_B21*TY_G111*TY_G28 - TY_B25*TY_G17*TY_G28 - TY_B24*TY_G18*TY_G28 - TY_B23*TY_G19*TY_G28 + 2*TY_B12*TY_G210*TY_G28 - TY_B21*TY_G110*TY_G29 - TY_B25*TY_G16*TY_G29 - TY_B24*TY_G17*TY_G29 - 
779        TY_B23*TY_G18*TY_G29 - TY_B22*TY_G19*TY_G29 + 2*TY_B14*TY_G27*TY_G29 + TY_B34*(2*(TY_G113*TY_G13 + TY_G112*TY_G14 + TY_G111*TY_G15 + TY_G110*TY_G16 + TY_G17*TY_G19) + pow(TY_G18,2)) + 
780        TY_B32*(2*(TY_G113*TY_G15 + TY_G112*TY_G16 + TY_G111*TY_G17 + TY_G110*TY_G18) + pow(TY_G19,2)) + TY_B14*pow(TY_G28,2) + TY_B12*pow(TY_G29,2);
781       
782        TY_w[13] = 2*TY_B32*(TY_G113*TY_G16 + TY_G112*TY_G17 + TY_G111*TY_G18 + TY_G110*TY_G19) + 2*TY_B34*(TY_G113*TY_G14 + TY_G112*TY_G15 + TY_G111*TY_G16 + TY_G110*TY_G17 + TY_G18*TY_G19) - TY_B21*TY_G110*TY_G210 - 
783        TY_B25*TY_G16*TY_G210 - TY_B24*TY_G17*TY_G210 - TY_B23*TY_G18*TY_G210 - TY_B22*TY_G19*TY_G210 - TY_B25*TY_G15*TY_G211 - TY_B24*TY_G16*TY_G211 - TY_B23*TY_G17*TY_G211 - TY_B22*TY_G18*TY_G211 - TY_B21*TY_G19*TY_G211 - 
784        TY_B25*TY_G14*TY_G212 - TY_B24*TY_G15*TY_G212 - TY_B23*TY_G16*TY_G212 - TY_B22*TY_G17*TY_G212 - TY_B21*TY_G18*TY_G212 - TY_B25*TY_G13*TY_G213 - TY_B24*TY_G14*TY_G213 - TY_B23*TY_G15*TY_G213 - TY_B22*TY_G16*TY_G213 - 
785        TY_B21*TY_G17*TY_G213 - TY_B25*TY_G113*TY_G23 - TY_B25*TY_G112*TY_G24 - TY_B24*TY_G113*TY_G24 + 2*TY_B14*TY_G213*TY_G24 - TY_B25*TY_G111*TY_G25 - TY_B24*TY_G112*TY_G25 - TY_B23*TY_G113*TY_G25 + 
786        2*TY_B14*TY_G212*TY_G25 - TY_G214*(TY_B24*TY_G13 + TY_B23*TY_G14 + TY_B22*TY_G15 + TY_B21*TY_G16 - 2*TY_B14*TY_G23 - 2*TY_B12*TY_G25) - TY_B25*TY_G110*TY_G26 - TY_B24*TY_G111*TY_G26 - TY_B23*TY_G112*TY_G26 - 
787        TY_B22*TY_G113*TY_G26 + 2*TY_B14*TY_G211*TY_G26 + 2*TY_B12*TY_G213*TY_G26 - TY_B24*TY_G110*TY_G27 - TY_B23*TY_G111*TY_G27 - TY_B22*TY_G112*TY_G27 - TY_B21*TY_G113*TY_G27 - TY_B25*TY_G19*TY_G27 + 
788        2*TY_B14*TY_G210*TY_G27 + 2*TY_B12*TY_G212*TY_G27 - TY_B23*TY_G110*TY_G28 - TY_B22*TY_G111*TY_G28 - TY_B21*TY_G112*TY_G28 - TY_B25*TY_G18*TY_G28 - TY_B24*TY_G19*TY_G28 + 2*TY_B12*TY_G211*TY_G28 - 
789        TY_B22*TY_G110*TY_G29 - TY_B21*TY_G111*TY_G29 - TY_B25*TY_G17*TY_G29 - TY_B24*TY_G18*TY_G29 - TY_B23*TY_G19*TY_G29 + 2*TY_B12*TY_G210*TY_G29 + 2*TY_B14*TY_G28*TY_G29;
790       
791        TY_w[14] = -(TY_B22*TY_G110*TY_G210) - TY_B21*TY_G111*TY_G210 - TY_B25*TY_G17*TY_G210 - TY_B24*TY_G18*TY_G210 - TY_B23*TY_G19*TY_G210 - TY_B21*TY_G110*TY_G211 - TY_B25*TY_G16*TY_G211 - TY_B24*TY_G17*TY_G211 - 
792        TY_B23*TY_G18*TY_G211 - TY_B22*TY_G19*TY_G211 - TY_B25*TY_G15*TY_G212 - TY_B24*TY_G16*TY_G212 - TY_B23*TY_G17*TY_G212 - TY_B22*TY_G18*TY_G212 - TY_B21*TY_G19*TY_G212 - TY_B25*TY_G14*TY_G213 - TY_B24*TY_G15*TY_G213 - 
793        TY_B23*TY_G16*TY_G213 - TY_B22*TY_G17*TY_G213 - TY_B21*TY_G18*TY_G213 - TY_B25*TY_G113*TY_G24 - TY_B25*TY_G112*TY_G25 - TY_B24*TY_G113*TY_G25 + 2*TY_B14*TY_G213*TY_G25 - TY_B25*TY_G111*TY_G26 - TY_B24*TY_G112*TY_G26 - 
794        TY_B23*TY_G113*TY_G26 + 2*TY_B14*TY_G212*TY_G26 - TY_G214*(TY_B25*TY_G13 + TY_B24*TY_G14 + TY_B23*TY_G15 + TY_B22*TY_G16 + TY_B21*TY_G17 - 2*TY_B14*TY_G24 - 2*TY_B12*TY_G26) - TY_B25*TY_G110*TY_G27 - 
795        TY_B24*TY_G111*TY_G27 - TY_B23*TY_G112*TY_G27 - TY_B22*TY_G113*TY_G27 + 2*TY_B14*TY_G211*TY_G27 + 2*TY_B12*TY_G213*TY_G27 - TY_B24*TY_G110*TY_G28 - TY_B23*TY_G111*TY_G28 - TY_B22*TY_G112*TY_G28 - 
796        TY_B21*TY_G113*TY_G28 - TY_B25*TY_G19*TY_G28 + 2*TY_B14*TY_G210*TY_G28 + 2*TY_B12*TY_G212*TY_G28 - TY_B23*TY_G110*TY_G29 - TY_B22*TY_G111*TY_G29 - TY_B21*TY_G112*TY_G29 - TY_B25*TY_G18*TY_G29 - TY_B24*TY_G19*TY_G29 + 
797        2*TY_B12*TY_G211*TY_G29 + TY_B32*(2*(TY_G113*TY_G17 + TY_G112*TY_G18 + TY_G111*TY_G19) + pow(TY_G110,2)) + 
798        TY_B34*(2*(TY_G113*TY_G15 + TY_G112*TY_G16 + TY_G111*TY_G17 + TY_G110*TY_G18) + pow(TY_G19,2)) + TY_B12*pow(TY_G210,2) + TY_B14*pow(TY_G29,2); 
799       
800        TY_w[15] = 2*TY_B34*(TY_G113*TY_G16 + TY_G112*TY_G17 + TY_G111*TY_G18 + TY_G110*TY_G19) + 2*TY_B32*(TY_G110*TY_G111 + TY_G113*TY_G18 + TY_G112*TY_G19) - TY_B23*TY_G110*TY_G210 - TY_B22*TY_G111*TY_G210 - 
801        TY_B21*TY_G112*TY_G210 - TY_B25*TY_G18*TY_G210 - TY_B24*TY_G19*TY_G210 - TY_B22*TY_G110*TY_G211 - TY_B21*TY_G111*TY_G211 - TY_B25*TY_G17*TY_G211 - TY_B24*TY_G18*TY_G211 - TY_B23*TY_G19*TY_G211 + 
802        2*TY_B12*TY_G210*TY_G211 - TY_B21*TY_G110*TY_G212 - TY_B25*TY_G16*TY_G212 - TY_B24*TY_G17*TY_G212 - TY_B23*TY_G18*TY_G212 - TY_B22*TY_G19*TY_G212 - TY_B25*TY_G15*TY_G213 - TY_B24*TY_G16*TY_G213 - 
803        TY_B23*TY_G17*TY_G213 - TY_B22*TY_G18*TY_G213 - TY_B21*TY_G19*TY_G213 - TY_B25*TY_G113*TY_G25 - TY_B25*TY_G112*TY_G26 - TY_B24*TY_G113*TY_G26 + 2*TY_B14*TY_G213*TY_G26 - TY_B25*TY_G111*TY_G27 - TY_B24*TY_G112*TY_G27 - 
804        TY_B23*TY_G113*TY_G27 + 2*TY_B14*TY_G212*TY_G27 - TY_G214*(TY_B25*TY_G14 + TY_B24*TY_G15 + TY_B23*TY_G16 + TY_B22*TY_G17 + TY_B21*TY_G18 - 2*TY_B14*TY_G25 - 2*TY_B12*TY_G27) - TY_B25*TY_G110*TY_G28 - 
805        TY_B24*TY_G111*TY_G28 - TY_B23*TY_G112*TY_G28 - TY_B22*TY_G113*TY_G28 + 2*TY_B14*TY_G211*TY_G28 + 2*TY_B12*TY_G213*TY_G28 - TY_B24*TY_G110*TY_G29 - TY_B23*TY_G111*TY_G29 - TY_B22*TY_G112*TY_G29 - 
806        TY_B21*TY_G113*TY_G29 - TY_B25*TY_G19*TY_G29 + 2*TY_B14*TY_G210*TY_G29 + 2*TY_B12*TY_G212*TY_G29;
807       
808        TY_w[16] = -(TY_B24*TY_G110*TY_G210) - TY_B23*TY_G111*TY_G210 - TY_B22*TY_G112*TY_G210 - TY_B21*TY_G113*TY_G210 - TY_B25*TY_G19*TY_G210 - TY_B23*TY_G110*TY_G211 - TY_B22*TY_G111*TY_G211 - TY_B21*TY_G112*TY_G211 - 
809        TY_B25*TY_G18*TY_G211 - TY_B24*TY_G19*TY_G211 - TY_B22*TY_G110*TY_G212 - TY_B21*TY_G111*TY_G212 - TY_B25*TY_G17*TY_G212 - TY_B24*TY_G18*TY_G212 - TY_B23*TY_G19*TY_G212 + 2*TY_B12*TY_G210*TY_G212 - 
810        TY_B21*TY_G110*TY_G213 - TY_B25*TY_G16*TY_G213 - TY_B24*TY_G17*TY_G213 - TY_B23*TY_G18*TY_G213 - TY_B22*TY_G19*TY_G213 - TY_B25*TY_G113*TY_G26 - TY_B25*TY_G112*TY_G27 - TY_B24*TY_G113*TY_G27 + 
811        2*TY_B14*TY_G213*TY_G27 - TY_B25*TY_G111*TY_G28 - TY_B24*TY_G112*TY_G28 - TY_B23*TY_G113*TY_G28 + 2*TY_B14*TY_G212*TY_G28 - 
812        TY_G214*(TY_B25*TY_G15 + TY_B24*TY_G16 + TY_B23*TY_G17 + TY_B22*TY_G18 + TY_B21*TY_G19 - 2*TY_B14*TY_G26 - 2*TY_B12*TY_G28) - TY_B25*TY_G110*TY_G29 - TY_B24*TY_G111*TY_G29 - TY_B23*TY_G112*TY_G29 - 
813        TY_B22*TY_G113*TY_G29 + 2*TY_B14*TY_G211*TY_G29 + 2*TY_B12*TY_G213*TY_G29 + TY_B34*(2*(TY_G113*TY_G17 + TY_G112*TY_G18 + TY_G111*TY_G19) + pow(TY_G110,2)) + 
814        TY_B32*(2*TY_G110*TY_G112 + 2*TY_G113*TY_G19 + pow(TY_G111,2)) + TY_B14*pow(TY_G210,2) + TY_B12*pow(TY_G211,2);
815       
816        TY_w[17] = 2*TY_B32*(TY_G111*TY_G112 + TY_G110*TY_G113) + 2*TY_B34*(TY_G110*TY_G111 + TY_G113*TY_G18 + TY_G112*TY_G19) - TY_B25*TY_G110*TY_G210 - TY_B24*TY_G111*TY_G210 - TY_B23*TY_G112*TY_G210 - TY_B22*TY_G113*TY_G210 - 
817        TY_B24*TY_G110*TY_G211 - TY_B23*TY_G111*TY_G211 - TY_B22*TY_G112*TY_G211 - TY_B21*TY_G113*TY_G211 - TY_B25*TY_G19*TY_G211 + 2*TY_B14*TY_G210*TY_G211 - TY_B23*TY_G110*TY_G212 - TY_B22*TY_G111*TY_G212 - 
818        TY_B21*TY_G112*TY_G212 - TY_B25*TY_G18*TY_G212 - TY_B24*TY_G19*TY_G212 + 2*TY_B12*TY_G211*TY_G212 - TY_B22*TY_G110*TY_G213 - TY_B21*TY_G111*TY_G213 - TY_B25*TY_G17*TY_G213 - TY_B24*TY_G18*TY_G213 - 
819        TY_B23*TY_G19*TY_G213 + 2*TY_B12*TY_G210*TY_G213 - TY_B25*TY_G113*TY_G27 - TY_B25*TY_G112*TY_G28 - TY_B24*TY_G113*TY_G28 + 2*TY_B14*TY_G213*TY_G28 - TY_B25*TY_G111*TY_G29 - TY_B24*TY_G112*TY_G29 - 
820        TY_B23*TY_G113*TY_G29 + 2*TY_B14*TY_G212*TY_G29 - TY_G214*(TY_B21*TY_G110 + TY_B25*TY_G16 + TY_B24*TY_G17 + TY_B23*TY_G18 + TY_B22*TY_G19 - 2*TY_B14*TY_G27 - 2*TY_B12*TY_G29);
821       
822        TY_w[18] = -(TY_B25*TY_G111*TY_G210) - TY_B24*TY_G112*TY_G210 - TY_B23*TY_G113*TY_G210 - TY_B25*TY_G110*TY_G211 - TY_B24*TY_G111*TY_G211 - TY_B23*TY_G112*TY_G211 - TY_B22*TY_G113*TY_G211 - TY_B24*TY_G110*TY_G212 - 
823        TY_B23*TY_G111*TY_G212 - TY_B22*TY_G112*TY_G212 - TY_B21*TY_G113*TY_G212 - TY_B25*TY_G19*TY_G212 + 2*TY_B14*TY_G210*TY_G212 - TY_B23*TY_G110*TY_G213 - TY_B22*TY_G111*TY_G213 - TY_B21*TY_G112*TY_G213 - 
824        TY_B25*TY_G18*TY_G213 - TY_B24*TY_G19*TY_G213 + 2*TY_B12*TY_G211*TY_G213 - TY_B25*TY_G113*TY_G28 - 
825        TY_G214*(TY_B22*TY_G110 + TY_B21*TY_G111 + TY_B25*TY_G17 + TY_B24*TY_G18 + TY_B23*TY_G19 - 2*TY_B12*TY_G210 - 2*TY_B14*TY_G28) - TY_B25*TY_G112*TY_G29 - TY_B24*TY_G113*TY_G29 + 2*TY_B14*TY_G213*TY_G29 + 
826        TY_B34*(2*TY_G110*TY_G112 + 2*TY_G113*TY_G19 + pow(TY_G111,2)) + TY_B32*(2*TY_G111*TY_G113 + pow(TY_G112,2)) + TY_B14*pow(TY_G211,2) + TY_B12*pow(TY_G212,2);
827       
828        TY_w[19] = 2*TY_B32*TY_G112*TY_G113 + 2*TY_B34*(TY_G111*TY_G112 + TY_G110*TY_G113) - TY_B25*TY_G112*TY_G210 - TY_B24*TY_G113*TY_G210 - TY_B25*TY_G111*TY_G211 - TY_B24*TY_G112*TY_G211 - TY_B23*TY_G113*TY_G211 - 
829        TY_B25*TY_G110*TY_G212 - TY_B24*TY_G111*TY_G212 - TY_B23*TY_G112*TY_G212 - TY_B22*TY_G113*TY_G212 + 2*TY_B14*TY_G211*TY_G212 - TY_B24*TY_G110*TY_G213 - TY_B23*TY_G111*TY_G213 - TY_B22*TY_G112*TY_G213 - 
830        TY_B21*TY_G113*TY_G213 - TY_B25*TY_G19*TY_G213 + 2*TY_B14*TY_G210*TY_G213 + 2*TY_B12*TY_G212*TY_G213 - TY_B25*TY_G113*TY_G29 - 
831        TY_G214*(TY_B23*TY_G110 + TY_B22*TY_G111 + TY_B21*TY_G112 + TY_B25*TY_G18 + TY_B24*TY_G19 - 2*TY_B12*TY_G211 - 2*TY_B14*TY_G29);
832       
833        TY_w[20] = -(TY_B25*TY_G113*TY_G210) - TY_B25*TY_G112*TY_G211 - TY_B24*TY_G113*TY_G211 - TY_B25*TY_G111*TY_G212 - TY_B24*TY_G112*TY_G212 - TY_B23*TY_G113*TY_G212 - TY_B25*TY_G110*TY_G213 - TY_B24*TY_G111*TY_G213 - 
834        TY_B23*TY_G112*TY_G213 - TY_B22*TY_G113*TY_G213 + 2*TY_B14*TY_G211*TY_G213 - (TY_B24*TY_G110 + TY_B23*TY_G111 + TY_B22*TY_G112 + TY_B21*TY_G113 + TY_B25*TY_G19 - 2*TY_B14*TY_G210 - 2*TY_B12*TY_G212)*TY_G214 + 
835        TY_B34*(2*TY_G111*TY_G113 + pow(TY_G112,2)) + TY_B32*pow(TY_G113,2) + TY_B14*pow(TY_G212,2) + TY_B12*pow(TY_G213,2);
836       
837        TY_w[21] = TY_B25*(TY_A23*TY_B14*(-3*TY_A52*TY_B24*TY_B25 + (2*TY_A43*TY_B24 + TY_A42*TY_B25)*TY_B34) + TY_B25*(TY_A22*TY_B14*(-(TY_A52*TY_B25) + TY_A43*TY_B34) 
838                                                                                                                                                         + TY_A12*(4*TY_A52*TY_B24*TY_B25 - (3*TY_A43*TY_B24 + TY_A42*TY_B25)*TY_B34)))*pow(TY_B34,3);
839       
840        TY_w[22] = (-(TY_A23*TY_B14) + TY_A12*TY_B25)*(TY_A52*TY_B25 - TY_A43*TY_B34)*pow(TY_B25,2)*pow(TY_B34,3);
841/*     
842        if( debug )
843        {
844                sprintf(buf,"\rCoefficients of polynomial\r");
845                XOPNotice(buf);
846                int i;
847                for ( i = 0; i < 23; i ++ ) {
848                        sprintf(buf, "w[%d] = %f\r", i, TY_w[i] );
849                        XOPNotice(buf);
850                }
851                sprintf(buf, "\r" );
852                XOPNotice(buf);
853        }
854 */
855}
856
857double TY_Q( double d2 )
858{
859        return d2 * TY_B32 + pow( d2, 3 ) *  TY_B34;
860}
861
862double TY_V( double d2 )
863{
864        return  -( pow( d2, 2 ) * TY_G13 + pow( d2, 3 ) * TY_G14 + pow( d2, 4 ) * TY_G15 + pow( d2, 5 ) * TY_G16 + 
865                          pow( d2, 6 ) * TY_G17 + pow( d2, 7 ) *  TY_G18 + pow( d2, 8 ) * TY_G19 + pow( d2, 9 ) * TY_G110 + 
866                          pow( d2, 10 ) *  TY_G111 + pow( d2, 11 ) *  TY_G112 + pow( d2, 12 ) *  TY_G113 );
867}
868
869double TY_W( double d2 )
870{
871        return d2  * TY_G22 + pow( d2, 2 ) * TY_G23 + pow( d2, 3 ) * TY_G24 + pow( d2, 4 ) * TY_G25 + pow( d2, 5 ) * TY_G26 + 
872        pow( d2, 6 ) * TY_G27 + pow( d2, 7 ) *  TY_G28 + pow( d2, 8 ) *  TY_G29 + pow( d2, 9 ) * TY_G210 + 
873        pow( d2, 10 ) * TY_G211 + pow( d2, 11 ) * TY_G212 + pow( d2, 12 ) * TY_G213 + pow( d2, 13 ) * TY_G214;
874}
875
876double TY_X( double d2 )
877{
878        return TY_V( d2 ) / TY_W( d2 );
879}
880
881// solve the linear system depending on d1, d2 using Cramer's rule
882void TY_SolveLinearEquations( double d1, double d2, 
883                                                      double* a, double* b, double* c1, double* c2)
884{
885        double det    = TY_q22 * d1 * d2;
886        double det_a  = TY_qa12  * d2 + TY_qa21  * d1 + TY_qa22  * d1 * d2 + TY_qa23  * d1 * pow( d2, 2 ) + TY_qa32  * pow( d1, 2 ) * d2; 
887        double det_b  = TY_qb12  * d2 + TY_qb21  * d1 + TY_qb22  * d1 * d2 + TY_qb23  * d1 * pow( d2, 2 ) + TY_qb32  * pow( d1, 2 ) * d2; 
888        double det_c1 = TY_qc112 * d2 + TY_qc121 * d1 + TY_qc122 * d1 * d2 + TY_qc123 * d1 * pow( d2, 2 ) + TY_qc132 * pow( d1, 2 ) * d2; 
889        double det_c2 = TY_qc212 * d2 + TY_qc221 * d1 + TY_qc222 * d1 * d2 + TY_qc223 * d1 * pow( d2, 2 ) + TY_qc232 * pow( d1, 2 ) * d2;
890       
891        *a  = det_a  / det;
892        *b  = det_b  / det;
893        *c1 = det_c1 / det;
894        *c2 = det_c2 / det;
895}
896
897//Solve the system of linear and nonlinear equations for given Zi, Ki, phi which gives at
898// most 22 solutions for the parameters a,b,ci,di. From the set of solutions choose the
899// physical one and return it.
900int TY_SolveEquations( double Z1, double Z2, double K1, double K2, double phi, 
901                                           double* a, double* b, double* c1, double* c2, double* d1, double* d2, 
902                                       int debug )
903{
904       
905        // the two coupled non-linear eqautions were reduced to a
906        // 22nd order polynomial, the roots are give all possible solutions
907        // for d2, than d1 can be computed by the function X
908        char buf[256];
909       
910        double real_coefficient[23];
911        double imag_coefficient[23];
912       
913        double real_root[22];
914        double imag_root[22];
915       
916        //integer degree of polynomial
917        int degree = 22;
918        int i;
919        double x,y;
920        double var_a, var_b, var_c1, var_c2, var_d1, var_d2;
921        double sol_a[22], sol_b[22], sol_c1[22], sol_c2[22], sol_d1[22], sol_d2[22];
922       
923        int j = 0;
924       
925        int n_roots,n,selected_root;
926        double dr,qmax,q,dq,min,sum;
927        double *sq,*gr;
928        double Pi = 3.14159265358979323846264338327950288;   /* pi */
929       
930       
931        // reduce system to a polynomial from which all solution are extracted
932        // by doing that a lot of global background variables are set
933        TY_ReduceNonlinearSystem( Z1, Z2, K1, K2, phi, debug );
934       
935        // vector of real and imaginary coefficients in order of decreasing powers
936        for ( i = 0; i <= degree; i++ )
937        {
938                // the global variablw TY_w was set by TY_ReduceNonlinearSystem
939                real_coefficient[i] = TY_w[ 22 - i ];
940                imag_coefficient[i] = 0.;
941        }
942       
943        // Jenkins-Traub method to approximate the roots of the polynomial
944        cpoly( real_coefficient, imag_coefficient, degree, real_root, imag_root );
945       
946        // show the result if in debug mode
947/*     
948        if ( debug )
949        {
950                for ( i = 0; i < degree; i++ )
951                {
952                        x = real_root[i];
953                        y = imag_root[i];
954                        if ( chop( y ) == 0 )
955                                sprintf(buf, "root(%d) = %g\r", i+1, x );
956                        else
957                                sprintf(buf, "root(%d) = %g + %g i\r", i+1, x, y );
958                        XOPNotice(buf);
959                }
960                sprintf(buf, "\r" );
961                XOPNotice(buf);
962        }
963 */
964       
965        // select real roots and those satisfying Q(x) != 0 and W(x) != 0
966        // Paper: Cluster formation in two-Yukawa Fluids, J. Chem. Phys. 122, 2005
967        // The right set of (a, b, c1, c2, d1, d2) should have the following properties:
968        // (1) a > 0
969        // (2) d1, d2 are real
970        // (3) vi/Ki > 0 <=> g(Zi) > 0
971        // (4) if there is still more than root, calculate g(r) for each root
972        //     and g(r) of the correct root should have the minimum average value
973        //         inside the hardcore 
974
975        for ( i = 0; i < degree; i++ )
976        {
977                x = real_root[i];
978                y = imag_root[i];
979               
980                if ( chop( y ) == 0 && TY_W( x ) != 0 && TY_Q( x ) != 0 )
981                {
982                        var_d1 = TY_X( x );
983                        var_d2 = x;
984                       
985                        // solution of linear system for given d1, d2 to obtain a,b,ci,di
986                        TY_SolveLinearEquations( var_d1, var_d2, &var_a, &var_b, &var_c1, &var_c2 );
987                       
988                        // select physical solutions, for details check paper: "Cluster formation in
989                        // two-Yukawa fluids", J. Chem. Phys. 122 (2005)
990                        if ( var_a > 0 && 
991                                 TY_g( Z1, phi, Z1, Z2, var_a, var_b, var_c1, var_c2, var_d1, var_d2 ) > 0 &&
992                                 TY_g( Z2, phi, Z1, Z2, var_a, var_b, var_c1, var_c2, var_d1, var_d2 ) > 0 )
993                        {
994                                sol_a[j]  = var_a;
995                                sol_b[j]  = var_b;
996                                sol_c1[j] = var_c1;
997                                sol_c2[j] = var_c2;
998                                sol_d1[j] = var_d1;
999                                sol_d2[j] = var_d2;
1000/*                             
1001                                if ( debug )
1002                                {
1003                                        double 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] ) );
1004                                        double 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] ) );
1005                                        double 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] ) );
1006                                        double 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] ) );
1007                                        double 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] ) );
1008                                        double 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] ) );
1009                                       
1010                                        sprintf(buf, "solution[%d] = (%g, %g, %g, %g, %g, %g), ( eq == 0 ) = (%g, %g, %g, %g, %g, %g)\r", j,
1011                                                   sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j],
1012                                                   eq1 , eq2, eq3, eq4, eq5, eq6 );
1013                                        XOPNotice(buf);
1014                                }
1015 */
1016                                j++;
1017                        }
1018                }
1019        }
1020        // number  remaining roots
1021        n_roots = j;
1022       
1023       
1024        // if there is still more than one root left, than choose the one with the minimum
1025        // average value inside the hardcore
1026        if ( n_roots > 1 )
1027        {
1028                // the number of q values should be a power of 2
1029                // in order to speed up the FFT
1030                n = 1 << 14;
1031               
1032                // the maximum q value should be large enough
1033                // to enable a reasoble approximation of g(r)
1034                qmax = 16 * 10 * 2 * Pi;
1035                dq = qmax / ( n - 1 );
1036               
1037                // step size for g(r) = dr
1038               
1039                // allocate memory for pair correlation function g(r)
1040                // and structure factor S(q)
1041                sq = malloc( sizeof( double ) * n );
1042                gr = malloc( sizeof( double ) * n );
1043               
1044                // loop over all remaining roots
1045                min = 1e99;
1046                selected_root = 10;     
1047                sum = 0;
1048                for ( j = 0; j < n_roots; j++) 
1049                {
1050                        // calculate structure factor at different q values
1051                        for ( i = 0; i < n; i++) 
1052                        {
1053                                q = dq * i;
1054                                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] );
1055/*                             
1056                                if(i<20 && debug) {
1057                                        sprintf(buf, "after SqTwoYukawa: s(q) = %g\r",sq[i] );
1058                                        XOPNotice(buf);
1059                                }
1060 */
1061                        }
1062                       
1063                        // calculate pair correlation function for given
1064                        // structure factor, g(r) is computed at values
1065                        // r(i) = i * dr
1066                        PairCorrelation( phi, dq, sq, &dr, gr, n );
1067                        // determine sum inside the hardcore
1068                        // 0 =< r < 1 of the pair-correlation function
1069                        sum = 0;
1070                        for (i = 0; i < floor( 1. / dr ); i++ )  {
1071                                sum += fabs( gr[i] );
1072/*                             
1073                                if(i<20 && debug) {
1074                                        sprintf(buf, "g(r) in core = %g\r",fabs(gr[i]));
1075                                        XOPNotice(buf);
1076                                }
1077*/                             
1078                        }
1079
1080                        if ( sum < min )
1081                        {
1082                                min = sum;
1083                                selected_root = j;
1084                        }
1085                }       
1086                free( gr );
1087                free( sq );
1088               
1089                // physical solution was found
1090                *a  = sol_a [ selected_root ];//sol_a [ selected_root ];
1091                *b  = sol_b [ selected_root ];
1092                *c1 = sol_c1[ selected_root ];
1093                *c2 = sol_c2[ selected_root ];
1094                *d1 = sol_d1[ selected_root ];
1095                *d2 = sol_d2[ selected_root ];
1096               
1097                return 1;
1098        }
1099        else if ( n_roots == 1 ) 
1100        {
1101                *a  = sol_a [0];
1102                *b  = sol_b [0];
1103                *c1 = sol_c1[0];
1104                *c2 = sol_c2[0];
1105                *d1 = sol_d1[0];
1106                *d2 = sol_d2[0];
1107               
1108                return 1;
1109        }
1110
1111        // no solution was found
1112        return 0;
1113}
Note: See TracBrowser for help on using the repository browser.