Index: sansmodels/src/libigor/2Y_OneYukawa.c
===================================================================
--- sansmodels/src/libigor/2Y_OneYukawa.c	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
+++ sansmodels/src/libigor/2Y_OneYukawa.c	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
@@ -0,0 +1,370 @@
+/*
+ *  Yukawa.c
+ *  twoyukawa
+ *
+ *  Created by Marcus Hennig on 5/12/10.
+ *  Copyright 2010 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "2Y_OneYukawa.h"
+#include "2Y_cpoly.h"
+#include "2Y_utility.h"
+#include "2Y_PairCorrelation.h"
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+
+/*
+================================================================================================== 
+ 
+ Structure factor for the potential 
+ 
+ V(r) = -kB * T * K * exp[ -Z * (r - 1)] / r for r > 1
+ V(r) = inf for r <= 1
+ 
+ The structure factor is parametrized by (a, b, c, d) 
+ which depend on (Z, K, phi).  
+ 
+================================================================================================== 
+*/
+
+double Y_sigma( double s, double Z, double a, double b, double c, double d )
+{
+	return -(a / 2. + b + c * exp( -Z )) / s + a * pow( s, -3 ) + b * pow( s, -2 ) + ( c + d ) * pow( s + Z, -1 );
+}
+
+double Y_tau( double s, double Z, double a, double b, double c )
+{
+	return b * pow( s, -2 ) + a * ( pow( s, -3 ) + pow( s, -2 ) ) - pow( s, -1 ) * c * Z * exp( -Z ) * pow( s + Z, -1 );
+} 
+
+double Y_q( double s, double Z, double a, double b, double c, double d )
+{
+	return Y_sigma(s, Z, a, b, c, d ) - exp( -s ) * Y_tau( s, Z, a, b, c );
+}
+
+double Y_g( double s, double phi, double Z, double a, double b, double c, double d )
+{
+	return s * Y_tau( s, Z, a, b, c ) * exp( -s ) / ( 1 - 12 * phi * Y_q( s, Z, a, b, c, d ) );
+}
+
+double Y_hq( double q, double Z, double K, double v )
+{
+	double t1, t2, t3, t4;
+
+	if ( q == 0) 
+	{
+		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.;
+	}
+	else 
+	{
+		
+		t1 = ( 1 - v / ( 2 * K * Z * exp( Z ) ) ) * ( ( 1 - cos( q ) ) / ( q*q ) - 1 / ( Z*Z + q*q ) );
+		t2 = ( v*v * ( q * cos( q ) - Z * sin( q ) ) ) / ( 4 * K * Z*Z * q * ( Z*Z + q*q ) );
+		t3 = ( q * cos( q ) + Z * sin( q ) ) / ( q * ( Z*Z + q*q ) );
+		t4 = v / ( Z * exp( Z ) ) - v*v / ( 4 * K * Z*Z * exp( 2 * Z ) ) - K;
+		
+		return v / Z * t1 - t2 + t3 * t4;
+	}
+}
+
+double Y_pc( double q, 
+			double Z, double K, double phi,
+			double a, double b, double c, double d )
+{	
+	double v = 24 * phi * K * exp( Z ) * Y_g( Z, phi, Z, a, b, c, d );
+	
+	double a0 = a * a;
+	double b0 = -12 * phi *( pow( a + b,2 ) / 2 + a * c * exp( -Z ) );
+	
+	double t1, t2, t3, t4;
+	
+	if ( q == 0 ) 
+	{
+		t1 = a0 / 3;
+		t2 = b0 / 4;
+		t3 = a0 * phi / 12;
+	}
+	else 
+	{
+		t1 = a0 * ( sin( q ) - q * cos( q ) ) / pow( q, 3 );
+		t2 = b0 * ( 2 * q * sin( q ) - ( q * q - 2 ) * cos( q ) - 2 ) / pow( q, 4 );
+		t3 = a0 * phi * ( ( q*q - 6 ) * 4 * q * sin( q ) - ( pow( q, 4 ) - 12 * q*q + 24) * cos( q ) + 24 ) / ( 2 * pow( q, 6 ) );
+	}
+	t4 = Y_hq( q, Z, K, v );
+	return -24 * phi * ( t1 + t2 + t3 + t4 );
+}
+
+double SqOneYukawa( double q, 
+				 double Z, double K, double phi,
+				 double a, double b, double c, double d )
+{
+	//structure factor one-yukawa potential
+	return 1. / ( 1. - Y_pc( q, Z, K, phi, a, b, c, d ) );
+}
+
+
+/*
+================================================================================================== 
+
+ The structure factor S(q) for the one-Yukawa potential is parameterized by a,b,c,d which are
+ the solution of a system of non-linear equations given Z,K, phi. There are at most 4 solutions 
+ from which the physical one is chosen
+ 
+================================================================================================== 
+ */
+
+double Y_LinearEquation_1( double Z, double K, double phi, double a, double b, double c, double d )
+{
+	return b - 12*phi*(-a/8. - b/6. + d*pow(Z,-2) + c*(pow(Z,-2) - exp(-Z)*(0.5 + (1 + Z)*pow(Z,-2))));
+}
+
+double Y_LinearEquation_2( double Z, double K, double phi, double a, double b, double c, double d )
+{
+	return 1 - a - 12*phi*(-a/3. - b/2. + d*pow(Z,-1) + c*(pow(Z,-1) - (1 + Z)*exp(-Z)*pow(Z,-1)));
+}
+
+double Y_LinearEquation_3( double Z, double K, double phi, double a, double b, double c, double d )
+{
+	return K*exp(Z) - Z*d*(1-12*phi*Y_q(Z, Z, a, b, c, d));
+}
+
+double Y_NonlinearEquation( double Z, double K, double phi, double a, double b, double c, double d )
+{
+	return c + d - 12*phi*((c + d)*Y_sigma(Z, Z, a, b, c, d) - c*exp(-Z)*Y_tau(Z, Z, a, b, c));
+}
+
+// Check the computed solutions satisfy the system of equations 
+int Y_CheckSolution( double Z, double K, double phi, 
+					 double a, double b, double c, double d )
+{
+	double eq_1 = chop( Y_LinearEquation_1 ( Z, K, phi, a, b, c, d ) );
+	double eq_2 = chop( Y_LinearEquation_2 ( Z, K, phi, a, b, c, d ) );
+	double eq_3 = chop( Y_LinearEquation_3 ( Z, K, phi, a, b, c, d ) );
+	double eq_4 = chop( Y_NonlinearEquation( Z, K,  phi, a, b, c, d ) );
+	
+	// check if all equation are zero
+	return eq_1 == 0 && eq_2 == 0 && eq_3 == 0 && eq_4 == 0;
+}
+
+int Y_SolveEquations( double Z, double K, double phi, double* a, double* b, double* c, double* d, int debug )
+{
+	char buf[256];
+
+	// at most there are 4 solutions for a,b,c,d
+	double sol_a[4], sol_b[4], sol_c[4], sol_d[4];
+	
+	double m11 = (3*phi)/2.;
+	double m13 = 6*phi*exp(-Z)*(2 + Z*(2 + Z) - 2*exp(Z))*pow(Z,-2);
+	double m23 = -12*phi*exp(-Z)*(-1 - Z + exp(Z))*pow(Z,-1);
+	double m31 = -6*phi*exp(-Z)*pow(Z,-2)*(2*(1 + Z) + exp(Z)*(-2 + pow(Z,2)));
+	double m32 = -12*phi*(-1 + Z + exp(-Z))*pow(Z,-1);
+	double m33 = 6*phi*exp(-2*Z)*pow(-1 + exp(Z),2);
+	
+	double delta = m23*m31 - m13*m32 + m11*(-4*m13*m31 + (4*m23*m31)/3. + (8*m13*m32)/3. - m23*m32) + m33 + (4*(-3 + m11)*m11*m33)/9.;
+	double a1 = -(K*(m23 + (4*m11*(-3*m13 + m23))/3.)*exp(Z));
+	double a2 = -(m13*(m32 + 4*m11*Z)) + ((3 + 4*m11)*(m33 + m23*Z))/3.;
+	double a3 = -2*phi*pow(Z,-2)*(6*m23*m32 - 24*m11*m33 + 2*Z*((3 + 4*m11)*m33 - 3*m13*(m32 + 2*m11*Z)) + (3 + 4*m11)*m23*pow(Z,2));
+	
+	double b1 = -(K*((-3 + 8*m11)*m13 - 3*m11*m23)*exp(Z))/3.;
+	double b2 = m13*(m31 - Z + (8*m11*Z)/3.) - m11*(m33 + m23*Z);
+	double b3 = 2*phi*pow(Z,-2)*(m13*Z*(-6*m31 + 3*Z - 8*m11*Z) + 2*m33*(3 - 8*m11 + 3*m11*Z) + 3*m23*(2*m31 + m11*pow(Z,2)));
+	
+	double c1 = -(K*exp(Z)*pow(3 - 2*m11,2))/9.;
+	double c2 = -((3 + 4*m11)*m31)/3. + m11*m32 + Z + (4*(-3 + m11)*m11*Z)/9.;
+	double c3 = (-2*phi*pow(Z,-2)*(6*(12*m11*m31 + 3*m32 - 8*m11*m32) - 6*((3 + 4*m11)*m31 - 3*m11*m32)*Z + pow(3 - 2*m11,2)*pow(Z,2)))/3.;
+	
+	// determine d, as roots of the polynomial, from that build a,b,c
+	
+	double real_coefficient[5];
+	double imag_coefficient[5];
+	
+	double real_root[4];
+	double imag_root[4];
+	
+	double zeta = 24*phi*pow(-6*phi*Z*cosh(Z/2.) + (12*phi + (-1 + phi)*pow(Z,2))*sinh(Z/2.),2);	
+	double A[5];
+	int degree,i,j,n_roots;
+	double x,y;
+	int n,selected_root;
+	double qmax,q,dq,min,sum,dr;
+	double *sq,*gr;
+	
+	
+	A[0] = -(exp(3*Z)*pow(K,2)*pow(-1 + phi,2)*pow(Z,3) / zeta );
+	A[1] = K*Z*exp(Z)*(6*phi*(2 + 4*phi + (2 + phi)*Z) + exp(Z)*
+					   ((-24 + Z*(18 + (-6 + Z)*Z))*pow(phi,2) - 2*phi*(6 + (-3 + Z)*pow(Z,2)) + pow(Z,3))) / zeta;
+	A[2] = -12*K*phi*exp(Z)*pow(-1 + phi,2)*pow(Z,3)/zeta;
+	A[3] = 6*phi*Z*exp(-Z)*(-12*phi*(1 + 2*phi)*(-1 + exp(Z)) + 6*phi*Z*(3*phi + (2 + phi)*exp(Z)) + 
+							6*(-1 + phi)*phi*pow(Z,2) + pow(-1 + phi,2)*pow(Z,3))/zeta;
+	A[4] = -36*exp(-Z)*pow(-1 + phi,2)*pow(phi,2)*pow(Z,3)/zeta;
+/*	
+	if ( debug )
+	{
+		sprintf (buf, "(Z,K,phi) = (%g, %g, %g)\r", K, Z, phi );
+		XOPNotice(buf);
+		sprintf (buf, "A = (%g, %g, %g, %g, %g)\r", A[0], A[1], A[2], A[3], A[4] );
+		XOPNotice(buf);
+	}
+*/	
+	//integer degree of polynomial
+	degree = 4;
+	
+	// vector of real and imaginary coefficients in order of decreasing powers
+	for ( i = 0; i <= degree; i++ )
+	{
+		real_coefficient[i] = A[4-i];
+		imag_coefficient[i] = 0.;
+	}
+	
+	// Jenkins-Traub method to approximate the roots of the polynomial
+	cpoly( real_coefficient, imag_coefficient, degree, real_root, imag_root );
+	
+	// show the result if in debug mode
+/*
+	if ( debug )
+	{
+		for ( i = 0; i < degree; i++ )
+		{
+			x = real_root[i];
+			y = imag_root[i];
+			sprintf(buf, "root(%d) = %g + %g i\r", i+1, x, y );
+			XOPNotice(buf);
+		}
+		sprintf(buf, "\r" );
+		XOPNotice(buf);
+	}
+ */
+	// determine the set of solutions for a,b,c,d,
+	j = 0;
+	for ( i = 0; i < degree; i++ ) 
+	{
+		x = real_root[i];
+		y = imag_root[i];
+		
+		if ( chop( y ) == 0 ) 
+		{
+			sol_a[j] = ( a1 + a2 * x + a3 * x * x ) / ( delta * x );
+			sol_b[j] = ( b1 + b2 * x + b3 * x * x ) / ( delta * x );
+			sol_c[j] = ( c1 + c2 * x + c3 * x * x ) / ( delta * x );
+			sol_d[j] = x;
+	
+			j++;
+		}
+	}
+
+	// number  remaining roots 
+	n_roots = j;
+	
+//	sprintf(buf, "inside Y_solveEquations OK, before malloc: n_roots = %d\r",n_roots);
+//	XOPNotice(buf);	
+	
+	// if there is still more than one root left, than choose the one with the minimum
+	// average value inside the hardcore
+	
+	
+	
+	if ( n_roots > 1 )
+	{
+		// the number of q values should be a power of 2
+		// in order to speed up the FFT
+		n = 1 << 14;		//= 16384
+		
+		// the maximum q value should be large enough 
+		// to enable a reasoble approximation of g(r)
+		qmax = 1000.;
+		dq = qmax / ( n - 1 );
+		
+		// step size for g(r) = dr
+		
+		// allocate memory for pair correlation function g(r)
+		// and structure factor S(q)
+		// (note that sq and gr are pointers)
+		sq = malloc( sizeof( double ) * n );
+		gr = malloc( sizeof( double ) * n );
+		
+		// loop over all remaining roots
+		min = 1e99;
+		selected_root=0;	
+		
+//		sprintf(buf, "after malloc: n,dq = %d  %g\r",n,dq);
+//		XOPNotice(buf);
+				
+		for ( j = 0; j < n_roots; j++) 
+		{
+			// calculate structure factor at different q values
+			for ( i = 0; i < n ; i++) 
+			{
+				q = dq * i;
+				sq[i] = SqOneYukawa( q, Z, K, phi, sol_a[j], sol_b[j], sol_c[j], sol_d[j] );
+/*				
+				if(i<20 && debug) {
+					sprintf(buf, "after SqOneYukawa: s(q) = %g\r",sq[i] );
+					XOPNotice(buf);	
+				}
+ */
+			}
+						
+			// calculate pair correlation function for given
+			// structure factor, g(r) is computed at values
+			// r(i) = i * dr
+			PairCorrelation( phi, dq, sq, &dr, gr, n );
+			
+//			sprintf(buf, "after PairCorrelation: \r");
+//			XOPNotice(buf);
+			
+			// determine sum inside the hardcore 
+			// 0 =< r < 1 of the pair-correlation function
+			sum = 0;
+			for (i = 0; i < floor( 1. / dr ); i++ ) 
+			{
+				sum += fabs( gr[i] );
+/*				
+				if(i<20 && debug) {
+					sprintf(buf, "g(r) in core = %g\r",fabs(gr[i]));
+					XOPNotice(buf);
+				}
+ */
+			}
+			
+//			sprintf(buf, "after hard core: sum, min = %g %g\r",sum,min);
+//			XOPNotice(buf);
+			
+			if ( sum < min )
+			{
+				min = sum;
+				selected_root = j;
+			}
+			
+			
+		}	
+		free( gr );
+		free( sq );
+		
+//		sprintf(buf, "after free: selected root = %d\r",selected_root);
+//		XOPNotice(buf);
+		
+		// physical solution was found
+		*a = sol_a[ selected_root ];
+		*b = sol_b[ selected_root ];
+		*c = sol_c[ selected_root ];
+		*d = sol_d[ selected_root ];
+		
+//		sprintf(buf, "after solution found (ret 1): \r");
+//		XOPNotice(buf);
+		
+		return 1;
+	}
+	else if ( n_roots == 1 ) 
+	{
+		*a = sol_a[0];
+		*b = sol_b[0];
+		*c = sol_c[0];
+		*d = sol_d[0];
+		
+		return 1;
+	}
+	// no solution was found
+	return 0;
+}
Index: sansmodels/src/libigor/2Y_OneYukawa.h
===================================================================
--- sansmodels/src/libigor/2Y_OneYukawa.h	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
+++ sansmodels/src/libigor/2Y_OneYukawa.h	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
@@ -0,0 +1,14 @@
+/*
+ *  Yukawa.h
+ *  twoyukawa
+ *
+ *  Created by Marcus Hennig on 5/12/10.
+ *  Copyright 2010 __MyCompanyName__. All rights reserved.
+ *
+ */
+int Y_CheckSolution( double Z, double K, double phi, double a, double b, double c, double d );
+int Y_SolveEquations( double Z, double K, double phi, 
+					  double* a, double* b, double* c, double* d, int debug );
+double SqOneYukawa( double q, 
+				double Z, double K, double phi,
+				double a, double b, double c, double d );
Index: sansmodels/src/libigor/2Y_PairCorrelation.c
===================================================================
--- sansmodels/src/libigor/2Y_PairCorrelation.c	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
+++ sansmodels/src/libigor/2Y_PairCorrelation.c	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
@@ -0,0 +1,226 @@
+/*
+ *  PairCorrelation.c
+ *  twoyukawa
+ *
+ *  Created by Marcus Hennig on 5/9/10.
+ *  Copyright 2010 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "2Y_PairCorrelation.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+//#include <gsl/gsl_errno.h> 
+//#include <gsl/gsl_fft_real.h>
+
+/*
+=================================================================================================== 
+
+ Source: J.B.Hayter: A Program for the fast bi-directional transforms
+ between g(r) and S(Q), ILL internal scientific report, October 1979
+ 
+ The transformation between structure factor and pair correlation 
+ function is given by 
+ 
+ g(x) = 1 + 1 / (12*pi*phi*x) * int( [S(q)-q]*q*sin(q*x), { 0, inf } )
+ 
+ where phi is the volume frcation, x and q are dimensionless variables, 
+ scaled by the radius a of the particles:
+ 
+ r = x * a;
+ Q = q / a
+ 
+ Discretizing the integral leads to
+ 
+ x[k] = 2*pi*k / (N * dq)
+ g(x[k]) = 1 + N*dq^3 / (24*pi^2*phi*k) * Im{ sum(S[n]*exp(2*pi*i*n*k/N),{n,0,N-1}) }
+ 
+ where S[n] = n*(S(q[n])-1) with q[n]=n * dq
+
+=================================================================================================== 
+*/
+
+/*
+int PairCorrelation_GSL( double phi, double dq, double* Sq, double* dr, double* gr, int N )
+{
+	double* data = malloc( sizeof(double) * N );
+	int n;
+	for ( n = 0; n < N; n++ )
+		data[n] = n * ( Sq[n] - 1 );
+	
+	// data[k] -> sum( data[n] * exp(-2*pi*i*n*k/N), {n, 0, N-1 })
+	int stride = 1;
+	int error  = gsl_fft_real_radix2_transform( data, stride, N );
+	
+	// if no errors detected 
+	if ( error == GSL_SUCCESS ) 
+	{
+		double alpha = N * pow( dq, 3 ) / ( 24 * M_PI * M_PI * phi );
+	
+	
+		*dr = 2 * M_PI / ( N * dq );  
+		int k;
+		double real, imag;
+		for ( k = 0; k < N; k++ )
+		{
+			// the solutions of the transform is stored in data, 
+			// consult GSL manual for more details
+			if ( k == 0 || k == N / 2)
+			{ 
+				real = data[k];
+				imag = 0;
+			}
+			else if ( k < N / 2 ) 
+			{
+				real = data[k];
+				imag = data[N-k];
+			}
+			else if ( k > N / 2 )
+			{
+				real =  data[N-k];
+				imag = -data[k];	
+			}
+			
+			if ( k == 0 ) 
+				gr[k] = 0;
+			else
+				gr[k] = 1. + alpha / k * (-imag);
+		} 
+	} 
+	// if N is not a power of two
+	else if ( error == GSL_EDOM )
+	{
+		printf( "N is not a power of 2\n" );
+	}
+	else
+	{
+		printf( "Could not perform DFT (discrete fourier transform)\n" );
+	}
+	// release allocated memory
+	free( data );
+	
+	// return error value
+	return error;
+}
+ */
+
+
+// this uses numerical recipes for the FFT
+//
+int PairCorrelation( double phi, double dq, double* Sq, double* dr, double* gr, int N )
+{
+	double* data = malloc( sizeof(double) * N * 2);
+	int n,error,k;
+	double alpha,real,imag;
+	double Pi = 3.14159265358979323846264338327950288;   /* pi */
+
+	
+	for ( n = 0; n < N; n++ ) {
+		data[2*n] = n * ( Sq[n] - 1 );
+		data[2*n+1] = 0;
+	}
+	//	printf("start of new fft\n");
+	
+	// data[k] -> sum( data[n] * exp(-2*pi*i*n*k/N), {n, 0, N-1 })
+	//	int error  = gsl_fft_real_radix2_transform( data, stride, N );
+	error  = 1;
+	dfour1( data-1, N, 1 );		//N is the number of complex points
+	
+	//	printf("dfour1 is done\n");
+	
+	// if no errors detected 
+	if ( error == 1 ) 
+	{
+		alpha = N * pow( dq, 3 ) / ( 24 * Pi * Pi * phi );
+		
+		*dr = 2 * Pi / ( N * dq );  
+		for ( k = 0; k < N; k++ )
+		{
+			// the solutions of the transform is stored in data, 
+			// consult GSL manual for more details
+			if ( 2*k == 0 || 2*k == 2*N / 2)
+			{ 
+				real = data[2*k];
+				imag = 0;
+			}
+			else if ( 2*k < 2*N / 2 ) 
+			{
+				real = data[2*k];
+				imag = data[2*k+1];
+			}
+			else if ( 2*k > 2*N / 2 )
+			{
+				real =  data[2*k];
+				imag = -data[2*k+1];	
+			}
+			
+			if ( k == 0 ) 
+				gr[k] = 0;
+			else
+//				gr[k] = 1. + alpha / k * (-imag);		//if using GSL
+				gr[k] = 1. + alpha / k * (imag);		//if using NR
+		} 
+	} 
+	
+	// release allocated memory
+	free( data );
+	
+//	printf(" done with FFT assignment -- Using Numerical Recipes, not GSL\n");
+	
+	// return error value
+	return error;
+}
+
+
+// isign == 1 means no scaling of output. isign == -1 multiplies output by nn
+//
+//
+#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
+
+void dfour1(double data[], unsigned long nn, int isign)
+{
+	unsigned long n,mmax,m,j,istep,i;
+	double wtemp,wr,wpr,wpi,wi,theta;
+	double tempr,tempi;
+	
+	n=nn << 1;
+	j=1;
+	for (i=1;i<n;i+=2) {
+		if (j > i) {
+			SWAP(data[j],data[i]);
+			SWAP(data[j+1],data[i+1]);
+		}
+		m=n >> 1;
+		while (m >= 2 && j > m) {
+			j -= m;
+			m >>= 1;
+		}
+		j += m;
+	}
+	mmax=2;
+	while (n > mmax) {
+		istep=mmax << 1;
+		theta=isign*(6.28318530717959/mmax);
+		wtemp=sin(0.5*theta);
+		wpr = -2.0*wtemp*wtemp;
+		wpi=sin(theta);
+		wr=1.0;
+		wi=0.0;
+		for (m=1;m<mmax;m+=2) {
+			for (i=m;i<=n;i+=istep) {
+				j=i+mmax;
+				tempr=wr*data[j]-wi*data[j+1];
+				tempi=wr*data[j+1]+wi*data[j];
+				data[j]=data[i]-tempr;
+				data[j+1]=data[i+1]-tempi;
+				data[i] += tempr;
+				data[i+1] += tempi;
+			}
+			wr=(wtemp=wr)*wpr-wi*wpi+wr;
+			wi=wi*wpr+wtemp*wpi+wi;
+		}
+		mmax=istep;
+	}
+}
+#undef SWAP
Index: sansmodels/src/libigor/2Y_PairCorrelation.h
===================================================================
--- sansmodels/src/libigor/2Y_PairCorrelation.h	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
+++ sansmodels/src/libigor/2Y_PairCorrelation.h	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
@@ -0,0 +1,13 @@
+/*
+ *  PairCorrelation.h
+ *  twoyukawa
+ *
+ *  Created by Marcus Hennig on 5/9/10.
+ *  Copyright 2010 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+//int PairCorrelation_GSL( double phi, double dq, double* Sq, double* dr, double* gr, int N );
+int PairCorrelation( double phi, double dq, double* Sq, double* dr, double* gr, int N );
+
+void dfour1(double data[], unsigned long nn, int isign);
Index: sansmodels/src/libigor/2Y_TwoYukawa.c
===================================================================
--- sansmodels/src/libigor/2Y_TwoYukawa.c	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
+++ sansmodels/src/libigor/2Y_TwoYukawa.c	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
@@ -0,0 +1,1113 @@
+/*
+ *  twoyukawa.c
+ *  twoyukawa
+ *
+ *  Created by Marcus Hennig on 5/7/10.
+ *  Copyright 2010 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "2Y_TwoYukawa.h"
+#include "2Y_cpoly.h"
+#include "2Y_utility.h"
+#include "2Y_PairCorrelation.h"
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+
+/*
+ ================================================================================================== 
+ 
+ The two-yukawa structure factor is uniquley determined by 6 parameters a, b, c1, c2, d1, d2,
+ which are the solution of a system of 6 equations ( 4 linear, 2 nonlinear ). The solution can 
+ constructed by the roots of a polynomial of 22nd degree. For more details see attached 
+ Mathematica snotebook, where a derivation is given
+ 
+ ================================================================================================== 
+ */
+
+double TY_q22;
+double TY_qa12, TY_qa21, TY_qa22, TY_qa23, TY_qa32;
+double TY_qb12, TY_qb21, TY_qb22, TY_qb23, TY_qb32;
+double TY_qc112, TY_qc121, TY_qc122, TY_qc123, TY_qc132;
+double TY_qc212, TY_qc221, TY_qc222, TY_qc223, TY_qc232;
+double TY_A12, TY_A21, TY_A22, TY_A23, TY_A32, TY_A41, TY_A42, TY_A43, TY_A52;
+double TY_B12, TY_B14, TY_B21, TY_B22, TY_B23, TY_B24, TY_B25, TY_B32, TY_B34;
+double 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;
+double 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;
+double TY_w[23];
+
+double TY_sigma( double s, 
+				 double Z1, double Z2, 
+			     double a, double b, double c1, double c2, double d1, double d2 )
+{
+	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 );
+}
+
+double TY_tau( double s, 
+			   double Z1, double Z2, 
+		       double a, double b, double c1, double c2 )
+{
+	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 ) );
+} 
+
+double TY_q( double s, 
+			     double Z1, double Z2, 
+			     double a, double b, double c1, double c2, double d1, double d2 )
+{
+	return TY_sigma(s, Z1, Z2, a, b, c1, c2, d1, d2) - exp( -s ) * TY_tau(s, Z1, Z2, a,b, c1, c2);
+}
+
+double TY_g( double s, 
+		     double phi, double Z1, double Z2, 
+		     double a, double b, double c1, double c2, double d1, double d2 )
+{
+	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 ) );
+}
+
+/*
+ ================================================================================================== 
+ 
+ Structure factor for the potential 
+ 
+ V(r) = -kB * T * ( K1 * exp[ -Z1 * (r - 1)] / r + K2 * exp[ -Z2 * (r - 1)] / r ) for r > 1
+ V(r) = inf for r <= 1
+ 
+ The structure factor is parametrized by (a, b, c1, c2, d1, d2) 
+ which depend on (K1, K2, Z1, Z2, phi).  
+ 
+ ================================================================================================== 
+ */
+
+double TY_hq( double q, double Z, double K, double v )
+{
+	double t1, t2, t3, t4;
+
+	if ( q == 0) 
+	{
+		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.;
+	}
+	else 
+	{
+		
+		t1 = ( 1 - v / ( 2 * K * Z * exp( Z ) ) ) * ( ( 1 - cos( q ) ) / ( q*q ) - 1 / ( Z*Z + q*q ) );
+		t2 = ( v*v * ( q * cos( q ) - Z * sin( q ) ) ) / ( 4 * K * Z*Z * q * ( Z*Z + q*q ) );
+		t3 = ( q * cos( q ) + Z * sin( q ) ) / ( q * ( Z*Z + q*q ) );
+		t4 = v / ( Z * exp( Z ) ) - v*v / ( 4 * K * Z*Z * exp( 2 * Z ) ) - K;
+		
+		return v / Z * t1 - t2 + t3 * t4;
+	}
+}
+
+double TY_pc( double q, 
+			  double Z1, double Z2,double K1, double K2, double phi,
+			  double a, double b, double c1, double c2, double d1, double d2 )
+{	
+	double v1 = 24 * phi * K1 * exp( Z1 ) * TY_g( Z1, phi, Z1, Z2, a, b, c1, c2, d1, d2 );
+	double v2 = 24 * phi * K2 * exp( Z2 ) * TY_g( Z2, phi, Z1, Z2, a, b, c1, c2, d1, d2 );
+	
+	double a0 = a * a;
+	double b0 = -12 * phi *( pow( a + b,2 ) / 2 + a * ( c1 * exp( -Z1 ) + c2 * exp( -Z2 ) ) );
+	
+	double t1, t2, t3,t4;
+	
+	if ( q == 0 ) 
+	{
+		t1 = a0 / 3;
+		t2 = b0 / 4;
+		t3 = a0 * phi / 12;
+	}
+	else 
+	{
+		t1 = a0 * ( sin( q ) - q * cos( q ) ) / pow( q, 3 );
+		t2 = b0 * ( 2 * q * sin( q ) - ( q * q - 2 ) * cos( q ) - 2 ) / pow( q, 4 );
+		t3 = a0 * phi * ( ( q*q - 6 ) * 4 * q * sin( q ) - ( pow( q, 4 ) - 12 * q*q + 24) * cos( q ) + 24 ) / ( 2 * pow( q, 6 ) );
+	}
+	t4 = TY_hq( q, Z1, K1, v1 ) + TY_hq( q, Z2, K2, v2 );
+	return -24 * phi * ( t1 + t2 + t3 + t4 );
+}
+
+double SqTwoYukawa( double q, 
+				    double Z1, double Z2, double K1, double K2, double phi,
+				    double a, double b, double c1, double c2, double d1, double d2 )
+{
+	if ( Z1 == Z2 ) 
+	{
+		// one-yukawa potential
+		return 0;
+	} 
+	else 
+	{
+		// two-yukawa potential
+		return 1. / ( 1. - TY_pc( q, Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) );
+	}
+}
+
+/*
+================================================================================================== 
+
+ Non-linear eqaution system that determines the parameter for structure factor
+  
+================================================================================================== 
+*/
+
+double TY_LinearEquation_1( double Z1, double Z2, double K1, double K2, double phi, 
+						    double a, double b, double c1, double c2, double d1, double d2 )
+{
+	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 ) ) ) );
+}
+
+double TY_LinearEquation_2( double Z1, double Z2, double K1, double K2, double phi, 
+							double a, double b, double c1, double c2, double d1, double d2 )
+{
+	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 ) ) );
+}	
+
+double TY_LinearEquation_3( double Z1, double Z2, double K1, double K2, double phi, 
+						    double a, double b, double c1, double c2, double d1, double d2 )
+{								
+	return K1 * exp( Z1 ) - d1 * Z1 * ( 1 - 12 * phi * TY_q( Z1, Z1, Z2, a, b, c1, c2, d1, d2 ) ); 
+}
+
+double TY_LinearEquation_4( double Z1, double Z2, double K1, double K2, double phi, 
+						    double a, double b, double c1, double c2, double d1, double d2 )
+{
+	return K2 * exp( Z2 ) - d2 * Z2 * ( 1 - 12 * phi * TY_q( Z2, Z1, Z2, a, b, c1, c2, d1, d2 ) );
+}
+
+double TY_NonlinearEquation_1( double Z1, double Z2, double K1, double K2, double phi, 
+						       double a, double b, double c1, double c2, double d1, double d2 )
+{
+	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 ) );
+}
+
+double TY_NonlinearEquation_2( double Z1, double Z2, double K1, double K2, double phi, 
+						       double a, double b, double c1, double c2, double d1, double d2 )
+{
+	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 ) );
+}
+
+// Check the computed solutions satisfy the system of equations 
+int TY_CheckSolution( double Z1, double Z2, double K1, double K2, double phi, 
+				      double a, double b, double c1, double c2, double d1, double d2 )
+{
+	double eq_1 = chop( TY_LinearEquation_1( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) );
+	double eq_2 = chop( TY_LinearEquation_2( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) );
+	double eq_3 = chop( TY_LinearEquation_3( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) );
+	double eq_4 = chop( TY_LinearEquation_4( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) );
+	double eq_5 = chop( TY_NonlinearEquation_1( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) );
+	double eq_6 = chop( TY_NonlinearEquation_2( Z1, Z2, K1, K2, phi, a, b, c1, c2, d1, d2 ) );
+	
+	// check if all equation are zero
+	return eq_1 == 0 && eq_2 == 0 && eq_3 == 0 && eq_4 == 0 && eq_5 == 0 && eq_6 == 0;
+}
+
+void TY_ReduceNonlinearSystem( double Z1, double Z2, double K1, double K2, double phi, int debug )
+{
+	/* solution of the 4 linear equations depending on d1 and d2, the solution is polynomial
+	 in d1, d2. We represend the solution as determiants obtained by Cramer's rule
+	 which can be expressed by their coefficient matrices
+	 */
+	char buf[256];
+	
+	double m11 = (3*phi)/2.;
+	double m13 = 6*phi*exp(-Z1)*(2 + Z1*(2 + Z1) - 2*exp(Z1))*pow(Z1,-2);
+	double m14 = 6*phi*exp(-Z2)*(2 + Z2*(2 + Z2) - 2*exp(Z2))*pow(Z2,-2);
+	double m23 = -12*phi*exp(-Z1)*(-1 - Z1 + exp(Z1))*pow(Z1,-1);
+	double m24 = -12*phi*exp(-Z2)*(-1 - Z2 + exp(Z2))*pow(Z2,-1);
+	double m31 = -6*phi*exp(-Z1)*pow(Z1,-2)*(2*(1 + Z1) + exp(Z1)*(-2 + pow(Z1,2)));
+	double m32 = -12*phi*(-1 + Z1 + exp(-Z1))*pow(Z1,-1);
+	double m33 = 6*phi*exp(-2*Z1)*pow(-1 + exp(Z1),2);
+	double m34 = 12*phi*exp(-Z1 - Z2)*(Z2 - (Z1 + Z2)*exp(Z1) + Z1*exp(Z1 + Z2))*pow(Z1 + Z2,-1);
+	double m41 = -6*phi*exp(-Z2)*pow(Z2,-2)*(2*(1 + Z2) + exp(Z2)*(-2 + pow(Z2,2)));
+	double m42 = -12*phi*(-1 + Z2 + exp(-Z2))*pow(Z2,-1);
+	double m43 = 12*phi*exp(-Z1 - Z2)*(Z1 - (Z1 + Z2 - Z2*exp(Z1))*exp(Z2))*pow(Z1 + Z2,-1);	
+	double m44 = 6*phi*exp(-2*Z2)*pow(-1 + exp(Z2),2);
+	
+	/* determinant of the linear system expressed as coefficient matrix in d1, d2 */
+	
+	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.) + 
+		  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*
+		  (m33*(3*m41 + 4*m11*m41 - 3*m11*m42) + (-3*m31 - 4*m11*m31 + 3*m11*m32)*m43) + 3*m23*(-3*m34*m41 - 4*m11*m34*m41 + 
+		  3*m11*m34*m42 + 3*m31*m44 + 4*m11*m31*m44 - 3*m11*m32*m44) - (m34*m43 - m33*m44)*pow(3 - 2*m11,2))/9.;
+/*	
+	if( debug ) 
+	{
+		sprintf(buf,"\rDet = \r" );
+		XOPNotice(buf);
+		sprintf(buf, "%f\t%f\r%f\t%f\r", 0., 0., 0., TY_q22 );
+		XOPNotice(buf);
+	}
+*/
+	
+	/* Matrix representation of the determinant of the of the system where row refering to 
+	 the variable a is replaced by solution vector */
+	
+	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.;
+	
+	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.;
+	
+	TY_qa22 = m14*(-(m23*m42*Z1) + 4*m11*m43*Z1 - m33*(m42 + 4*m11*Z2) + m32*(m43 + m23*Z2)) + 
+		   (3*m13*(m24*m42*Z1 - 4*m11*m44*Z1 + m34*(m42 + 4*m11*Z2) - m32*(m44 + m24*Z2)) + 
+		   (3 + 4*m11)*(-(m24*m43*Z1) + m23*m44*Z1 - m34*(m43 + m23*Z2) + m33*(m44 + m24*Z2)))/3.;	
+	
+	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)) + 
+		   3*m32*(Z1 + Z2)*(2*m43 + m13*pow(Z2,2))) + 
+		   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)) - 
+		   3*m32*(Z1 + Z2)*(2*m44 + m14*pow(Z2,2))) + 
+		   2*(3*(m14*m33*m42 - m13*m34*m42 - m14*m32*m43 + m34*m43 + m13*m32*m44 - m33*m44)*Z2*(Z1 + Z2) + 
+		   2*m11*(6*(-(m14*m43) + m13*m44)*Z1*pow(Z2,2) + m34*(Z1 + Z2)*(2*m43*(-3 + Z2) - 3*m13*pow(Z2,2)) + 
+		  m33*(Z1 + Z2)*(6*m44 - 2*m44*Z2 + 3*m14*pow(Z2,2)))))*pow(Z1 + Z2,-1);	
+	
+	TY_qa32 = 2*phi*pow(Z1,-2)*(m24*((-3*m13*m42 + (3 + 4*m11)*m43)*(Z1 + Z2)*pow(Z1,2) - 
+		   2*m33*(3*m42*(Z1 + Z2) + (3 + 4*m11)*Z2*pow(Z1,2)) + 6*m32*(m43*(Z1 + Z2) + m13*Z2*pow(Z1,2))) + 
+		   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)) - 
+		   6*m32*(m44*(Z1 + Z2) + m14*Z2*pow(Z1,2))) + 
+		   2*(3*(m14*m33*m42 - m13*m34*m42 - m14*m32*m43 + m34*m43 + m13*m32*m44 - m33*m44)*Z1*(Z1 + Z2) + 
+		   2*m11*(-3*(m14*m43 - m13*m44)*(Z1 + Z2)*pow(Z1,2) + 2*m34*(m43*(-3 + Z1)*(Z1 + Z2) - 3*m13*Z2*pow(Z1,2)) + 
+		  m33*(-2*m44*(-3 + Z1)*(Z1 + Z2) + 6*m14*Z2*pow(Z1,2)))))*pow(Z1 + Z2,-1);
+/*		
+	if( debug ) 
+	{
+		sprintf(buf,"\rDet_a = \r" );
+		XOPNotice(buf);
+		sprintf(buf, "%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. );
+		XOPNotice(buf);
+	}
+*/
+	
+	/* Matrix representation of the determinant of the of the system where row refering to 
+	 the variable b is replaced by solution vector */
+	
+	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.;
+	
+	TY_qb21 = (K2*(-3*m13*m24*m31 + 3*m11*m24*m33 + m14*(3*m23*m31 + (3 - 8*m11)*m33) - 3*m13*m34 + 8*m11*m13*m34 - 3*m11*m23*m34)*
+			exp(Z2))/3.;
+	
+	TY_qb22 = m13*(m31*m44 - m24*m41*Z1 - m44*Z1 + (8*m11*m44*Z1)/3. + m24*m31*Z2 + m34*(-m41 + Z2 - (8*m11*Z2)/3.)) + 
+		   m14*(m23*m41*Z1 + m43*Z1 - (8*m11*m43*Z1)/3. + m33*(m41 - Z2 + (8*m11*Z2)/3.) - m31*(m43 + m23*Z2)) + 
+		   m11*(m24*m43*Z1 - m23*m44*Z1 + m34*(m43 + m23*Z2) - m33*(m44 + m24*Z2));	
+	
+	TY_qb23 = 2*phi*(3*m14*m23*m31 - 3*m13*m24*m31 + 3*m14*m33 - 8*m11*m14*m33 + 3*m11*m24*m33 - 3*m13*m34 + 8*m11*m13*m34 - 
+		   3*m11*m23*m34 + 2*(3*m24*(m33*m41 - m31*m43) + m23*(-3*m34*m41 + 3*m31*m44) + (-3 + 8*m11)*(m34*m43 - m33*m44))*
+		   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) + 
+		   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*
+		   pow(Z1 + Z2,-1));
+	
+	TY_qb32 = 2*phi*pow(Z1,-2)*(6*(-(m34*(m23*m41 + m43)) + m24*(m33*m41 - m31*m43) + (m23*m31 + m33)*m44) + 
+		   6*(-(m14*m33*m41) + m13*m34*m41 + m14*m31*m43 - m13*m31*m44)*Z1 + 
+		   3*(m14*(2*m23*m31 + 2*m33 - m23*m41 - m43) + m13*(-2*m34 + m24*(-2*m31 + m41) + m44))*pow(Z1,2) + 
+		   (m11*Z2*(16*m34*m43 - 16*m33*m44 - 6*m34*m43*Z1 + 6*m33*m44*Z1 + 
+		   (6*m24*m33 - 3*m24*m43 + 8*m14*(-2*m33 + m43) + (8*m13 - 3*m23)*(2*m34 - m44))*pow(Z1,2)) + 
+		   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)))*
+		   pow(Z1 + Z2,-1) + 6*(-(m14*(m23*m31 + m33)) + m13*(m24*m31 + m34))*pow(Z1,3)*pow(Z1 + Z2,-1));
+/*		
+	if( debug ) 
+	{
+		sprintf(buf,"\rDet_b = \r" );
+		XOPNotice(buf);
+		sprintf(buf, "%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. );
+		XOPNotice(buf);
+	}
+*/
+	
+	/* Matrix representation of the determinant of the of the system where row refering to 
+	 the variable c1 is replaced by solution vector */
+	
+	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.;
+	
+	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.;
+	
+	TY_qc122 = m14*(-4*m11*m41*Z1 - m42*Z1 + (8*m11*m42*Z1)/3. + m32*(-m41 + Z2 - (8*m11*Z2)/3.) + m31*(m42 + 4*m11*Z2)) + 
+			(3*m34*((3 + 4*m11)*m41 - 3*m11*m42) + 9*m11*m32*m44 + 9*m24*m41*Z1 + 12*m11*m24*m41*Z1 - 9*m11*m24*m42*Z1 + 9*m44*Z1 - 
+			12*m11*m44*Z1 + 9*m11*m24*m32*Z2 - 3*(3 + 4*m11)*m31*(m44 + m24*Z2) - m34*Z2*pow(3 - 2*m11,2) + 4*m44*Z1*pow(m11,2))/9.;
+	
+	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)) - 
+			2*(m14*m42 - m44)*Z2*(-(Z1*Z2) + m31*(Z1 + Z2))) + 4*(-2*m44*Z1 + m34*(Z1 + Z2))*pow(m11,2)*pow(Z2,2) - 
+			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)) - 
+			m31*(Z1 + Z2)*(6*m42 + (3 + 4*m11)*pow(Z2,2))) - 
+			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) + 
+			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)) - 
+			4*m14*m32*pow(Z2,3)))*pow(Z1 + Z2,-1))/3.;
+	
+	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))) + 
+			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) + 
+			3*m24*(((3 + 4*m11)*m41 - 3*m11*m42)*(Z1 + Z2)*pow(Z1,2) + 6*m32*(m41*(Z1 + Z2) + m11*Z2*pow(Z1,2)) - 
+			2*m31*(3*m42*(Z1 + Z2) + (3 + 4*m11)*Z2*pow(Z1,2))) + 
+			6*m11*(Z1*(-8*m32*m44 + m34*(m42*(8 - 3*Z1) + 4*m41*(-3 + Z1)) - 4*m31*m44*(-3 + Z1) + 3*m32*m44*Z1 - 
+			2*(3*m14*m41 - 2*m14*m42 + m44)*pow(Z1,2)) + 
+			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) + 
+			m34*(m42*(8 - 3*Z1) + 4*m41*(-3 + Z1) + 4*pow(Z1,2)))))*pow(Z1 + Z2,-1))/3.;
+/*		
+	if( debug ) 
+	{
+		sprintf(buf,"\rDet_c1 = \r" );
+		XOPNotice(buf);
+		sprintf(buf, "%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. );
+		XOPNotice(buf);
+	}
+ */
+	/* Matrix representation of the determinant of the of the system where row refering to 
+	 the variable c1 is replaced by solution vector */
+	
+	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.;
+	
+	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.;
+	
+	TY_qc222 = m13*(4*m11*m41*Z1 + m42*Z1 - (8*m11*m42*Z1)/3. + m32*(m41 - Z2 + (8*m11*Z2)/3.) - m31*(m42 + 4*m11*Z2)) + 
+			(9*m31*m43 - 9*(m23*m41 + m43)*Z1 + 9*m23*m31*Z2 + 3*m11*
+			((-4*m23*m41 + 3*m23*m42 + 4*m43)*Z1 + 4*m31*(m43 + m23*Z2) - 3*m32*(m43 + m23*Z2)) + 
+			m33*(-3*(3 + 4*m11)*m41 + 9*m11*m42 + Z2*pow(3 - 2*m11,2)) - 4*m43*Z1*pow(m11,2))/9.;
+	
+	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)) + 
+			2*(m13*m42 - m43)*Z2*(-(Z1*Z2) + m31*(Z1 + Z2))) - 4*(-2*m43*Z1 + m33*(Z1 + Z2))*pow(m11,2)*pow(Z2,2) + 
+			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)) - 
+			m31*(Z1 + Z2)*(6*m42 + (3 + 4*m11)*pow(Z2,2))) + 
+			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) + 
+			m33*(Z1 + Z2)*(8*m42 + 4*m41*(-3 + Z2) - 3*m42*Z2 + 2*pow(Z2,2)) + 2*m31*(Z1 + Z2)*(6*m43 - 2*m43*Z2 + 3*m13*pow(Z2,2)) - 
+			4*m13*m32*pow(Z2,3)))*pow(Z1 + Z2,-1))/3.;
+	
+	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))) + 
+			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) + 
+			3*m23*(((3 + 4*m11)*m41 - 3*m11*m42)*(Z1 + Z2)*pow(Z1,2) + 6*m32*(m41*(Z1 + Z2) + m11*Z2*pow(Z1,2)) - 
+			2*m31*(3*m42*(Z1 + Z2) + (3 + 4*m11)*Z2*pow(Z1,2))) + 
+			6*m11*(Z1*(-8*m32*m43 + m33*(m42*(8 - 3*Z1) + 4*m41*(-3 + Z1)) - 4*m31*m43*(-3 + Z1) + 3*m32*m43*Z1 - 
+			2*(3*m13*m41 - 2*m13*m42 + m43)*pow(Z1,2)) + 
+			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) + 
+			m33*(m42*(8 - 3*Z1) + 4*m41*(-3 + Z1) + 4*pow(Z1,2)))))*pow(Z1 + Z2,-1))/3.;
+/*		
+	if( debug ) 
+	{
+		sprintf(buf,"\rDet_c2 = \r" );
+		XOPNotice(buf);
+		sprintf(buf, "%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. );
+		XOPNotice(buf);
+	}	
+*/
+	
+	/* coefficient matrices of nonlinear equation 1 */
+	
+	TY_A12 = 6*phi*TY_qc112*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*TY_qc212*exp(Z1)*(-Z2 + (Z1 + Z2)*exp(Z1))*pow(Z1,2) + 
+		  exp(Z2)*(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))) - 
+		  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);
+	
+	TY_A21 = 6*phi*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*(TY_qc121*TY_qc212 + TY_qc112*TY_qc221)*exp(Z1)*(-Z2 + (Z1 + Z2)*exp(Z1))*pow(Z1,2) + 
+		  exp(Z2)*(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) + 
+		  exp(2*Z1)*(2*Z1*(TY_qb21*TY_qc112*(-1 + Z1)*(Z1 + Z2) + TY_qb12*TY_qc121*(-1 + Z1)*(Z1 + Z2) - 
+		  Z1*(TY_qc121*TY_qc212*Z1 + TY_qc112*(TY_qc121 + TY_qc221)*Z1 + TY_qc112*TY_qc121*Z2)) + TY_qa21*TY_qc112*(Z1 + Z2)*(-2 + pow(Z1,2)) + 
+		  TY_qa12*TY_qc121*(Z1 + Z2)*(-2 + pow(Z1,2))) - 2*TY_qc112*TY_qc121*(Z1 + Z2)*pow(Z1,2)))*pow(Z1 + Z2,-1);
+	
+	TY_A22 = exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(12*phi*(TY_qc122*TY_qc212 + TY_qc112*TY_qc222)*exp(Z1)*(-Z2 + (Z1 + Z2)*exp(Z1))*pow(Z1,2) + 
+		  exp(Z2)*(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) - 
+		  2*phi*TY_qc112*TY_qc122*(Z1 + Z2)*pow(Z1,2) + exp(2*Z1)*
+		  (6*phi*(2*Z1*(TY_qb22*TY_qc112*(-1 + Z1)*(Z1 + Z2) + TY_qb12*TY_qc122*(-1 + Z1)*(Z1 + Z2) - 
+		  Z1*(TY_qc122*TY_qc212*Z1 + TY_qc112*(TY_qc122 + TY_qc222)*Z1 + TY_qc112*TY_qc122*Z2)) + TY_qa22*TY_qc112*(Z1 + Z2)*(-2 + pow(Z1,2)) + 
+		  TY_qa12*TY_qc122*(Z1 + Z2)*(-2 + pow(Z1,2))) + TY_q22*TY_qc112*(Z1 + Z2)*pow(Z1,3))))*pow(Z1 + Z2,-1);
+	
+	TY_A23 = 6*phi*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*(TY_qc123*TY_qc212 + TY_qc112*TY_qc223)*exp(Z1)*(-Z2 + (Z1 + Z2)*exp(Z1))*pow(Z1,2) + 
+		  exp(Z2)*(2*(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) + 
+		  exp(2*Z1)*(2*Z1*(TY_qb23*TY_qc112*(-1 + Z1)*(Z1 + Z2) + TY_qb12*TY_qc123*(-1 + Z1)*(Z1 + Z2) - 
+		  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)) + 
+		  TY_qa12*TY_qc123*(Z1 + Z2)*(-2 + pow(Z1,2))) - 2*TY_qc112*TY_qc123*(Z1 + Z2)*pow(Z1,2)))*pow(Z1 + Z2,-1);
+	
+	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) + 
+		  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 + 
+		  2*(TY_qc121*TY_qc123 + TY_qc112*TY_qc132)*Z1 + TY_qc122*(TY_qa22 + TY_qb22 + TY_qc122*Z1)))*(Z1 + Z2)*exp(Z1 + Z2) - 
+		  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) + 
+		  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) - 
+		  6*phi*(Z1 + Z2)*exp(Z2)*(2*TY_qc121*TY_qc123 + 2*TY_qc112*TY_qc132 + pow(TY_qc122,2))*pow(Z1,2) + 
+		  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)) + 
+		  TY_qa12*(Z1 + Z2)*(-2 + pow(Z1,2))) + TY_qc122*(Z1 + Z2)*pow(Z1,3)) + 
+		  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)) + 
+		  TY_qa23*TY_qc121*(Z1 + Z2)*(-2 + pow(Z1,2)) - 2*TY_qb32*TY_qc112*pow(Z1,2) - 2*TY_qb23*TY_qc121*pow(Z1,2) - 2*TY_qb22*TY_qc122*pow(Z1,2) - 
+		  2*TY_qb21*TY_qc123*pow(Z1,2) - 2*TY_qb12*TY_qc132*pow(Z1,2) + 
+		  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 + 
+		  (2*TY_qb32*TY_qc112 + 2*TY_qb23*TY_qc121 + (TY_qa22 + 2*TY_qb22 - TY_qc122)*TY_qc122 + (TY_qa21 + 2*TY_qb21 - 2*TY_qc121)*TY_qc123 + 
+		  (TY_qa12 + 2*TY_qb12 - 2*TY_qc112)*TY_qc132)*pow(Z1,2)) + 2*TY_qb32*TY_qc112*pow(Z1,3) + 2*TY_qb23*TY_qc121*pow(Z1,3) + TY_qa22*TY_qc122*pow(Z1,3) + 
+		  2*TY_qb22*TY_qc122*pow(Z1,3) + TY_qa21*TY_qc123*pow(Z1,3) + 2*TY_qb21*TY_qc123*pow(Z1,3) - 2*TY_qc121*TY_qc123*pow(Z1,3) + TY_qa12*TY_qc132*pow(Z1,3) + 
+		  2*TY_qb12*TY_qc132*pow(Z1,3) - 2*TY_qc112*TY_qc132*pow(Z1,3) - 2*TY_qc132*TY_qc212*pow(Z1,3) - 2*TY_qc123*TY_qc221*pow(Z1,3) - 
+		  2*TY_qc122*TY_qc222*pow(Z1,3) - 2*TY_qc121*TY_qc223*pow(Z1,3) - 2*TY_qc112*TY_qc232*pow(Z1,3) - pow(TY_qc122,2)*pow(Z1,3))))*pow(Z1 + Z2,-1);
+	
+	TY_A41 = 6*phi*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*exp(Z1)*
+		  ((-(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) + 
+		  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)*
+		  exp(Z1) - 2*TY_qc121*TY_qc132*(Z1 + Z2)*pow(Z1,2) + 
+		  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 + 
+		  (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)) + 
+		  Z1*(-2*(TY_qa32*TY_qc121 + TY_qa21*(TY_q22 + TY_qc132)) - 2*(TY_qb32*TY_qc121 + TY_qb21*(TY_q22 + TY_qc132))*Z1 + 
+		  (TY_qa32*TY_qc121 + 2*TY_qb32*TY_qc121 + TY_qa21*TY_qc132 + 2*TY_qb21*TY_qc132 - 2*TY_qc121*TY_qc132 + TY_q22*(TY_qa21 + 2*TY_qb21 - 2*TY_qc121 - 2*TY_qc221) - 
+		  2*TY_qc132*TY_qc221 - 2*TY_qc121*TY_qc232)*pow(Z1,2)))))*pow(Z1 + Z2,-1);
+	
+	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) + 
+		  Z1*(TY_qb32*TY_qc122 + (TY_qa22 + TY_qb22)*TY_qc132 + TY_qc122*(TY_q22 + 2*TY_qc132)*Z1))*(Z1 + Z2)*exp(Z1 + Z2) - 
+		  12*phi*(TY_qc132*TY_qc222 + TY_qc122*TY_qc232)*Z2*exp(Z1)*pow(Z1,2) + 
+		  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) + 
+		  exp(2*Z1 + Z2)*(6*phi*(2*Z1*(TY_qb32*TY_qc122*(-1 + Z1)*(Z1 + Z2) + TY_qb22*TY_qc132*(-1 + Z1)*(Z1 + Z2) - 
+		  Z1*(TY_qc132*TY_qc222*Z1 + TY_qc122*(TY_qc132 + TY_qc232)*Z1 + TY_qc122*TY_qc132*Z2)) + TY_qa32*TY_qc122*(Z1 + Z2)*(-2 + pow(Z1,2)) + 
+		  TY_qa22*TY_qc132*(Z1 + Z2)*(-2 + pow(Z1,2))) + (Z1 + Z2)*pow(TY_q22,2)*pow(Z1,3) + 
+		  TY_q22*(6*phi*(2*Z1*(TY_qb22*(-1 + Z1)*(Z1 + Z2) - Z1*((TY_qc122 + TY_qc222)*Z1 + TY_qc122*Z2)) + TY_qa22*(Z1 + Z2)*(-2 + pow(Z1,2))) + 
+		  TY_qc132*(Z1 + Z2)*pow(Z1,3))))*pow(Z1 + Z2,-1);
+	
+	TY_A43 = -6*phi*exp(-2*Z1 - Z2)*pow(TY_q22,-2)*pow(Z1,-3)*(2*exp(Z1)*
+		  ((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) + 
+		  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)*
+		  exp(Z1) + 2*TY_qc123*TY_qc132*(Z1 + Z2)*pow(Z1,2) + 
+		  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 - 
+		  (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)) + 
+		  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 + 
+		  (-(TY_qa32*TY_qc123) - (TY_qa23 + 2*TY_qb23)*TY_qc132 + TY_q22*(-TY_qa23 + 2*(-TY_qb23 + TY_qc123 + TY_qc132 + TY_qc223)) + 
+		  2*(-(TY_qb32*TY_qc123) + TY_qc132*(TY_qc123 + TY_qc223) + TY_qc123*TY_qc232) + 2*pow(TY_q22,2))*pow(Z1,2)))))*pow(Z1 + Z2,-1);
+	
+	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) + 
+		  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)) - 
+		  TY_qa32*(Z1 + Z2)*(-2 + pow(Z1,2))) + (Z1 + Z2)*pow(TY_qc132,2)*pow(Z1,2) - 
+		  2*TY_qc132*(Z1 + Z2)*exp(Z1)*(TY_qa32 + (TY_qa32 + TY_qb32)*Z1 + (TY_q22 + TY_qc132)*pow(Z1,2))))*pow(Z1 + Z2,-1);
+	
+	// normalize A
+	/*double norm_A = sqrt(pow(TY_A52,2)+pow(TY_A43,2)+pow(TY_A42, 2)+pow(TY_A41, 2)+pow(TY_A32, 2)+
+						 pow(TY_A23,2)+pow(TY_A22,2)+pow(TY_A21, 2)+pow(TY_A12, 2));
+	TY_A12 /= norm_A;
+	TY_A21 /= norm_A;
+	TY_A22 /= norm_A;
+	TY_A23 /= norm_A;
+	TY_A32 /= norm_A;
+	TY_A41 /= norm_A;
+	TY_A42 /= norm_A;
+	TY_A43 /= norm_A;
+	TY_A52 /= norm_A;*/
+/*	
+	if( debug ) 
+	{
+		sprintf(buf,"\rNonlinear equation 1 = \r" );
+		XOPNotice(buf);
+		sprintf(buf, "%f\t\t%f\t\t%f\r", 0.,   TY_A12, 0.   );
+		XOPNotice(buf);
+		sprintf(buf, "%f\t\t%f\t\t%f\r", TY_A21, TY_A22, TY_A23 );
+		XOPNotice(buf);
+		sprintf(buf, "%f\t\t%f\t\t%f\r",  0.,  TY_A32, 0.   );
+		XOPNotice(buf);
+		sprintf(buf, "%f\t\t%f\t\t%f\r", TY_A41, TY_A42, TY_A43 );
+		XOPNotice(buf);
+		sprintf(buf, "%f\t\t%f\t\t%f\r", 0.,   TY_A52, 0. );		
+		XOPNotice(buf);
+	}
+ */
+	
+	/* coefficient matrices of nonlinear equation 2 */
+	
+	TY_B12 = 6*phi*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-2*TY_qc212*TY_qc221*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 
+		  2*exp(Z2)*((Z1 + Z2)*(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) + 
+		  (-(TY_qc121*TY_qc212) - TY_qc112*TY_qc221)*Z1*pow(Z2,2)) + 
+		  exp(2*Z2)*(exp(Z1)*(2*Z2*(TY_qb21*TY_qc212*(-1 + Z2)*(Z1 + Z2) + TY_qb12*TY_qc221*(-1 + Z2)*(Z1 + Z2) - 
+		  Z2*(TY_qc212*TY_qc221*Z1 + TY_qc112*TY_qc221*Z2 + TY_qc212*(TY_qc121 + TY_qc221)*Z2)) + TY_qa21*TY_qc212*(Z1 + Z2)*(-2 + pow(Z2,2)) + 
+		  TY_qa12*TY_qc221*(Z1 + Z2)*(-2 + pow(Z2,2))) + 2*(TY_qc121*TY_qc212 + TY_qc112*TY_qc221)*(Z1 + Z2)*pow(Z2,2)))*pow(Z1 + Z2,-1);
+	
+	TY_B14 = 6*phi*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-2*TY_qc212*TY_qc223*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 
+		  2*exp(Z2)*((Z1 + Z2)*(TY_qa12*TY_qc223 + TY_qa23*TY_qc212*(1 + Z2) + Z2*(TY_qb23*TY_qc212 + (TY_qa12 + TY_qb12)*TY_qc223 + TY_qc212*(TY_q22 + 2*TY_qc223)*Z2))*
+		  exp(Z1) + (-(TY_qc123*TY_qc212) - TY_qc112*TY_qc223)*Z1*pow(Z2,2)) + 
+		  exp(2*Z2)*(2*(TY_qc123*TY_qc212 + TY_qc112*(TY_q22 + TY_qc223))*(Z1 + Z2)*pow(Z2,2) + 
+		  exp(Z1)*(-2*(TY_qa23*TY_qc212 + TY_qa12*(TY_q22 + TY_qc223))*Z1 - 
+		  2*(TY_qa23*TY_qc212 + TY_qa12*(TY_q22 + TY_qc223) + (TY_qb23*TY_qc212 + TY_qb12*(TY_q22 + TY_qc223))*Z1)*Z2 + 
+		  (-2*(TY_qb23*TY_qc212 + TY_qb12*(TY_q22 + TY_qc223)) + (TY_q22*(TY_qa12 + 2*TY_qb12 - 2*TY_qc212) + TY_qc212*(TY_qa23 + 2*TY_qb23 - 2*TY_qc223) + 
+		  (TY_qa12 + 2*TY_qb12)*TY_qc223)*Z1)*pow(Z2,2) + 
+		  (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)*
+		  pow(Z2,3))))*pow(Z1 + Z2,-1);
+	
+	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)) + 
+		  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))) + 
+		  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)))
+		  )*pow(Z1 + Z2,-1);
+	
+	TY_B22 = exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-12*phi*TY_qc221*TY_qc222*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 
+		  12*phi*exp(Z2)*((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) + 
+		  (-(TY_qc122*TY_qc221) - TY_qc121*TY_qc222)*Z1*pow(Z2,2)) + 
+		  exp(2*Z2)*(12*phi*(TY_qc122*TY_qc221 + TY_qc121*TY_qc222)*(Z1 + Z2)*pow(Z2,2) + 
+		  exp(Z1)*(6*phi*(2*Z2*(TY_qb22*TY_qc221*(-1 + Z2)*(Z1 + Z2) + TY_qb21*TY_qc222*(-1 + Z2)*(Z1 + Z2) - 
+		  Z2*(TY_qc221*TY_qc222*Z1 + TY_qc121*TY_qc222*Z2 + TY_qc221*(TY_qc122 + TY_qc222)*Z2)) + TY_qa22*TY_qc221*(Z1 + Z2)*(-2 + pow(Z2,2)) + 
+		  TY_qa21*TY_qc222*(Z1 + Z2)*(-2 + pow(Z2,2))) + TY_q22*TY_qc221*(Z1 + Z2)*pow(Z2,3))))*pow(Z1 + Z2,-1);
+	
+	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) + 
+		  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) + 
+		  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 + 
+		  Z2*(TY_q22*TY_qc221 + 2*TY_qc221*TY_qc223 + 2*TY_qc212*TY_qc232 + pow(TY_qc222,2)))) + 
+		  (-(TY_qc132*TY_qc212) - TY_qc123*TY_qc221 - TY_qc122*TY_qc222 - TY_qc121*TY_qc223 - TY_qc112*TY_qc232)*Z1*pow(Z2,2)) + 
+		  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) + 
+		  exp(Z1)*(TY_q22*TY_qc222*(Z1 + Z2)*pow(Z2,3) - 6*phi*
+		  (2*(TY_qa32*TY_qc212 + TY_qa23*TY_qc221 + TY_qa22*TY_qc222 + TY_qa21*(TY_q22 + TY_qc223) + TY_qa12*TY_qc232)*Z1 + 
+		  2*(TY_qa32*TY_qc212 + TY_qa23*TY_qc221 + TY_qa22*TY_qc222 + TY_qa21*(TY_q22 + TY_qc223) + TY_qa12*TY_qc232 + 
+		  (TY_qb32*TY_qc212 + TY_qb23*TY_qc221 + TY_qb22*TY_qc222 + TY_qb21*(TY_q22 + TY_qc223) + TY_qb12*TY_qc232)*Z1)*Z2 - 
+		  (-2*(TY_qb32*TY_qc212 + TY_qb23*TY_qc221 + TY_qb22*TY_qc222 + TY_qb21*(TY_q22 + TY_qc223) + TY_qb12*TY_qc232) + 
+		  (TY_q22*(TY_qa21 + 2*TY_qb21 - 2*TY_qc221) + (TY_qa22 + 2*TY_qb22 - TY_qc222)*TY_qc222 + TY_qc221*(TY_qa23 + 2*TY_qb23 - 2*TY_qc223) + TY_qa21*TY_qc223 + 
+		  2*TY_qb21*TY_qc223 + TY_qc212*(TY_qa32 + 2*TY_qb32 - 2*TY_qc232) + TY_qa12*TY_qc232 + 2*TY_qb12*TY_qc232)*Z1)*pow(Z2,2) - 
+		  (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 + 
+		  TY_qc221*(TY_qa23 + 2*TY_qb23 - 2*TY_qc123 - 2*TY_qc223) + (TY_qa21 + 2*TY_qb21 - 2*TY_qc121)*TY_qc223 + 
+		  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);
+	
+	TY_B24 = exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(12*phi*(Z1 + Z2)*
+		  (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) - 
+		  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) + 
+		  12*phi*(TY_qc123*TY_qc222 + TY_qc122*(TY_q22 + TY_qc223))*(Z1 + Z2)*exp(2*Z2)*pow(Z2,2) + 
+		  exp(Z1 + 2*Z2)*(6*phi*(2*Z2*(TY_qb23*TY_qc222*(-1 + Z2)*(Z1 + Z2) + TY_qb22*TY_qc223*(-1 + Z2)*(Z1 + Z2) - 
+		  Z2*(TY_qc222*TY_qc223*Z1 + TY_qc122*TY_qc223*Z2 + TY_qc222*(TY_qc123 + TY_qc223)*Z2)) + TY_qa23*TY_qc222*(Z1 + Z2)*(-2 + pow(Z2,2)) + 
+		  TY_qa22*TY_qc223*(Z1 + Z2)*(-2 + pow(Z2,2))) + (Z1 + Z2)*pow(TY_q22,2)*pow(Z2,3) + 
+		  TY_q22*(6*phi*(2*Z2*(TY_qb22*(-1 + Z2)*(Z1 + Z2) - Z2*(TY_qc222*Z1 + (TY_qc122 + TY_qc222)*Z2)) + TY_qa22*(Z1 + Z2)*(-2 + pow(Z2,2))) + 
+		  TY_qc223*(Z1 + Z2)*pow(Z2,3))))*pow(Z1 + Z2,-1);
+	
+	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) + 
+		  (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)) - 
+		  TY_qa23*(Z1 + Z2)*(-2 + pow(Z2,2))) - 2*TY_qc123*(Z1 + Z2)*pow(Z2,2)) + 
+		  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);
+	
+	TY_B32 = 6*phi*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(-2*TY_qc221*TY_qc232*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 
+		  2*exp(Z2)*((Z1 + Z2)*(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) + 
+		  (-(TY_qc132*TY_qc221) - TY_qc121*TY_qc232)*Z1*pow(Z2,2)) + 
+		  exp(2*Z2)*(exp(Z1)*(2*Z2*(TY_qb32*TY_qc221*(-1 + Z2)*(Z1 + Z2) + TY_qb21*TY_qc232*(-1 + Z2)*(Z1 + Z2) - 
+		  Z2*(TY_qc221*TY_qc232*Z1 + TY_qc121*TY_qc232*Z2 + TY_qc221*(TY_q22 + TY_qc132 + TY_qc232)*Z2)) + TY_qa32*TY_qc221*(Z1 + Z2)*(-2 + pow(Z2,2)) + 
+		  TY_qa21*TY_qc232*(Z1 + Z2)*(-2 + pow(Z2,2))) + 2*(TY_qc132*TY_qc221 + TY_qc121*TY_qc232)*(Z1 + Z2)*pow(Z2,2)))*pow(Z1 + Z2,-1);
+	
+	TY_B34 = -6*phi*exp(-Z1 - 2*Z2)*pow(TY_q22,-2)*pow(Z2,-3)*(2*TY_qc223*TY_qc232*(Z1 + Z2)*exp(Z1)*pow(Z2,2) + 
+		  2*exp(Z2)*(-((Z1 + Z2)*(TY_qa23*TY_qc232 + TY_qa32*TY_qc223*(1 + Z2) + Z2*(TY_qb32*TY_qc223 + TY_qc232*(TY_qa23 + TY_qb23 + TY_q22*Z2 + 2*TY_qc223*Z2)))*exp(Z1)) + 
+		  (TY_qc132*TY_qc223 + TY_qc123*TY_qc232)*Z1*pow(Z2,2)) + exp(2*Z2)*
+		  (-2*(TY_qc132*(TY_q22 + TY_qc223) + TY_qc123*TY_qc232)*(Z1 + Z2)*pow(Z2,2) + 
+		  exp(Z1)*(2*(TY_qa32*(TY_q22 + TY_qc223) + TY_qa23*TY_qc232)*Z1 + 
+		  2*(TY_qa32*(TY_q22 + TY_qc223) + TY_qa23*TY_qc232 + (TY_qb32*(TY_q22 + TY_qc223) + TY_qb23*TY_qc232)*Z1)*Z2 - 
+		  (-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)*
+		  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);
+	
+	/*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));
+	
+	TY_B12 /= norm_B;
+	TY_B14 /= norm_B;
+	TY_B21 /= norm_B;
+	TY_B22 /= norm_B;
+	TY_B23 /= norm_B;
+	TY_B24 /= norm_B;
+	TY_B25 /= norm_B;
+	TY_B32 /= norm_B;
+	TY_B34 /= norm_B; */
+/*	
+	if( debug ) 
+	{
+		sprintf(buf,"\rNonlinear equation 2 = \r" );
+		XOPNotice(buf);
+		sprintf(buf, "%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0.,  TY_B12, 0.,  TY_B14, 0.  );
+		XOPNotice(buf);
+		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 );
+		XOPNotice(buf);
+		sprintf(buf, "%f\t\t%f\t\t%f\t\t%f\t\t%f\r", 0.,  TY_B32, 0.,  TY_B34, 0.  );
+		XOPNotice(buf);
+	}
+*/
+	
+	/* decrease order of nonlinear equation 1 by means of equation 2 */
+	
+	TY_F14 = -(TY_A32*TY_B12*TY_B32) + TY_A52*pow(TY_B12,2) + TY_A12*pow(TY_B32,2);
+	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;
+	TY_F18 = -(TY_A32*TY_B14*TY_B34) + TY_A52*pow(TY_B14,2) + TY_A12*pow(TY_B34,2);
+	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);
+	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);
+	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);
+	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;
+	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);
+	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);
+	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);
+	TY_F32 = -(TY_A41*TY_B21*TY_B32) + TY_A52*pow(TY_B21,2);
+	TY_F33 = 2*TY_A52*TY_B21*TY_B22 - TY_A42*TY_B21*TY_B32 - TY_A41*TY_B22*TY_B32;
+	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);
+	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;
+	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);
+	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;
+	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);
+	TY_F39 = 2*TY_A52*TY_B24*TY_B25 - TY_A43*TY_B24*TY_B34 - TY_A42*TY_B25*TY_B34;
+	TY_F310 = -(TY_A43*TY_B25*TY_B34) + TY_A52*pow(TY_B25,2);
+
+/*
+	if( debug ) 
+	{
+		sprintf(buf,"\rF = \r" );
+		XOPNotice(buf);
+		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.    );
+		XOPNotice(buf);
+		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.    );
+		XOPNotice(buf);
+		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 );
+		XOPNotice(buf);
+	}
+*/
+ 
+	TY_G13  = -(TY_B12*TY_F32);
+	TY_G14  = -(TY_B12*TY_F33);
+	TY_G15  = TY_B32*TY_F14 - TY_B14*TY_F32 - TY_B12*TY_F34;
+	TY_G16  = -(TY_B14*TY_F33) - TY_B12*TY_F35;
+	TY_G17  = TY_B34*TY_F14 + TY_B32*TY_F16 - TY_B14*TY_F34 - TY_B12*TY_F36;
+	TY_G18  = -(TY_B14*TY_F35) - TY_B12*TY_F37;
+	TY_G19  = TY_B34*TY_F16 + TY_B32*TY_F18 - TY_B14*TY_F36 - TY_B12*TY_F38;
+	TY_G110 = -(TY_B14*TY_F37) - TY_B12*TY_F39;
+	TY_G111 = TY_B34*TY_F18 - TY_B12*TY_F310 - TY_B14*TY_F38;
+	TY_G112 = -(TY_B14*TY_F39);
+	TY_G113 = -(TY_B14*TY_F310);
+	TY_G22  = -(TY_B21*TY_F32);
+	TY_G23  = -(TY_B22*TY_F32) - TY_B21*TY_F33;
+	TY_G24  = TY_B32*TY_F23 - TY_B23*TY_F32 - TY_B22*TY_F33 - TY_B21*TY_F34;
+	TY_G25  = TY_B32*TY_F24 - TY_B24*TY_F32 - TY_B23*TY_F33 - TY_B22*TY_F34 - TY_B21*TY_F35;
+	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;
+	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;
+	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;
+	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;
+	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;
+	TY_G211 = TY_B34*TY_F28 - TY_B22*TY_F310 - TY_B25*TY_F37 - TY_B24*TY_F38 - TY_B23*TY_F39;
+	TY_G212 = TY_B34*TY_F29 - TY_B23*TY_F310 - TY_B25*TY_F38 - TY_B24*TY_F39;
+	TY_G213 = -(TY_B24*TY_F310) - TY_B25*TY_F39;
+	TY_G214 = -(TY_B25*TY_F310);
+/*	
+	if( debug ) 
+	{
+		sprintf(buf,"\rG = \r" );
+		XOPNotice(buf);
+		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. );
+		XOPNotice(buf);
+		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 );
+		XOPNotice(buf);
+	}
+*/	
+	// coefficients for polynomial
+	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);
+	
+	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 - 
+	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;
+	
+	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 - 
+	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 + 
+	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);
+	
+	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 - 
+	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 - 
+	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 + 
+	2*TY_B12*TY_G22*TY_G27;
+	
+	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 - 
+	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 - 
+	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 + 
+	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) + 
+	TY_B12*pow(TY_G25,2);
+	
+	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 - 
+	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 - 
+	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 - 
+	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 - 
+	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;
+	
+	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 - 
+	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 - 
+	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 + 
+	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 - 
+	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)) + 
+	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);
+	
+	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 - 
+	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 - 
+	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 - 
+	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 + 
+	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 - 
+	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 + 
+	2*TY_B12*TY_G24*TY_G29;
+	
+	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 - 
+	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 + 
+	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 - 
+	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 - 
+	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 - 
+	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 + 
+	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) + 
+	TY_B12*pow(TY_G27,2);
+	
+	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 - 
+	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 - 
+	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 - 
+	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 + 
+	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 - 
+	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 - 
+	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 - 
+	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;
+	
+	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 - 
+	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 + 
+	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 - 
+	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 - 
+	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 - 
+	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 - 
+	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 + 
+	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)) + 
+	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);
+	
+	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 - 
+	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 - 
+	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 + 
+	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 - 
+	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 - 
+	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 + 
+	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 - 
+	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 + 
+	2*TY_B14*TY_G26*TY_G29 + 2*TY_B12*TY_G28*TY_G29;
+	
+	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 - 
+	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 - 
+	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 + 
+	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 - 
+	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 + 
+	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 - 
+	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 - 
+	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)) + 
+	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);
+	
+	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 - 
+	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 - 
+	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 - 
+	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 + 
+	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 - 
+	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 + 
+	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 - 
+	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;
+	
+	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 - 
+	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 - 
+	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 - 
+	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 - 
+	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 - 
+	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 + 
+	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)) + 
+	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); 
+	
+	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 - 
+	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 + 
+	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 - 
+	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 - 
+	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 - 
+	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 - 
+	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;
+	
+	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 - 
+	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 - 
+	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 + 
+	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 - 
+	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 - 
+	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)) + 
+	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);
+	
+	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 - 
+	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 - 
+	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 - 
+	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 - 
+	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);
+	
+	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 - 
+	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 - 
+	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 - 
+	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 + 
+	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);
+	
+	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 - 
+	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 - 
+	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 - 
+	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);
+	
+	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 - 
+	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 + 
+	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);
+	
+	TY_w[21] = TY_B25*(TY_A23*TY_B14*(-3*TY_A52*TY_B24*TY_B25 + (2*TY_A43*TY_B24 + TY_A42*TY_B25)*TY_B34) + TY_B25*(TY_A22*TY_B14*(-(TY_A52*TY_B25) + TY_A43*TY_B34) 
+																			 + TY_A12*(4*TY_A52*TY_B24*TY_B25 - (3*TY_A43*TY_B24 + TY_A42*TY_B25)*TY_B34)))*pow(TY_B34,3);
+	
+	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);
+/*	
+	if( debug ) 
+	{
+		sprintf(buf,"\rCoefficients of polynomial\r");
+		XOPNotice(buf);
+		int i;
+		for ( i = 0; i < 23; i ++ ) {
+			sprintf(buf, "w[%d] = %f\r", i, TY_w[i] );
+			XOPNotice(buf);
+		}
+		sprintf(buf, "\r" );
+		XOPNotice(buf);
+	}
+ */
+}
+
+double TY_Q( double d2 )
+{
+	return d2 * TY_B32 + pow( d2, 3 ) *  TY_B34;
+}
+
+double TY_V( double d2 )
+{
+	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 );
+}
+
+double TY_W( double d2 )
+{
+	return d2  * TY_G22 + pow( d2, 2 ) * TY_G23 + pow( d2, 3 ) * TY_G24 + pow( d2, 4 ) * TY_G25 + pow( d2, 5 ) * TY_G26 + 
+	pow( d2, 6 ) * TY_G27 + pow( d2, 7 ) *  TY_G28 + pow( d2, 8 ) *  TY_G29 + pow( d2, 9 ) * TY_G210 + 
+	pow( d2, 10 ) * TY_G211 + pow( d2, 11 ) * TY_G212 + pow( d2, 12 ) * TY_G213 + pow( d2, 13 ) * TY_G214;
+}
+
+double TY_X( double d2 )
+{
+	return TY_V( d2 ) / TY_W( d2 );
+}
+
+// solve the linear system depending on d1, d2 using Cramer's rule
+void TY_SolveLinearEquations( double d1, double d2, 
+						      double* a, double* b, double* c1, double* c2)
+{
+	double det    = TY_q22 * d1 * d2;
+	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; 
+	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; 
+	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; 
+	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;
+	
+	*a  = det_a  / det;
+	*b  = det_b  / det;
+	*c1 = det_c1 / det;
+	*c2 = det_c2 / det;
+}
+
+//Solve the system of linear and nonlinear equations for given Zi, Ki, phi which gives at 
+// most 22 solutions for the parameters a,b,ci,di. From the set of solutions choose the 
+// physical one and return it.
+int TY_SolveEquations( double Z1, double Z2, double K1, double K2, double phi, 
+					   double* a, double* b, double* c1, double* c2, double* d1, double* d2, 
+				       int debug )
+{
+	
+	// the two coupled non-linear eqautions were reduced to a
+	// 22nd order polynomial, the roots are give all possible solutions 
+	// for d2, than d1 can be computed by the function X 
+	char buf[256];
+	
+ 	double real_coefficient[23];
+	double imag_coefficient[23];
+	
+	double real_root[22];
+	double imag_root[22];
+	
+	//integer degree of polynomial
+	int degree = 22;
+	int i;
+	double x,y;
+	double var_a, var_b, var_c1, var_c2, var_d1, var_d2;
+	double sol_a[22], sol_b[22], sol_c1[22], sol_c2[22], sol_d1[22], sol_d2[22];
+	
+	int j = 0;
+	
+	int n_roots,n,selected_root;
+	double dr,qmax,q,dq,min,sum;
+	double *sq,*gr;
+	double Pi = 3.14159265358979323846264338327950288;   /* pi */
+	
+	
+	// reduce system to a polynomial from which all solution are extracted
+	// by doing that a lot of global background variables are set
+	TY_ReduceNonlinearSystem( Z1, Z2, K1, K2, phi, debug );
+	
+	// vector of real and imaginary coefficients in order of decreasing powers
+	for ( i = 0; i <= degree; i++ )
+	{
+		// the global variablw TY_w was set by TY_ReduceNonlinearSystem
+		real_coefficient[i] = TY_w[ 22 - i ];
+		imag_coefficient[i] = 0.;
+	}
+	
+	// Jenkins-Traub method to approximate the roots of the polynomial
+	cpoly( real_coefficient, imag_coefficient, degree, real_root, imag_root );
+	
+	// show the result if in debug mode
+/*	
+	if ( debug )
+	{
+		for ( i = 0; i < degree; i++ )
+		{
+			x = real_root[i];
+			y = imag_root[i];
+			if ( chop( y ) == 0 )
+				sprintf(buf, "root(%d) = %g\r", i+1, x );
+			else
+				sprintf(buf, "root(%d) = %g + %g i\r", i+1, x, y );
+			XOPNotice(buf);
+		}
+		sprintf(buf, "\r" );
+		XOPNotice(buf);
+	}
+ */
+	
+	// select real roots and those satisfying Q(x) != 0 and W(x) != 0
+	// Paper: Cluster formation in two-Yukawa Fluids, J. Chem. Phys. 122, 2005
+	// The right set of (a, b, c1, c2, d1, d2) should have the following properties:
+	// (1) a > 0 
+	// (2) d1, d2 are real
+	// (3) vi/Ki > 0 <=> g(Zi) > 0
+	// (4) if there is still more than root, calculate g(r) for each root
+	//     and g(r) of the correct root should have the minimum average value 
+	//	   inside the hardcore	
+
+	for ( i = 0; i < degree; i++ )
+	{
+		x = real_root[i];
+		y = imag_root[i];
+		
+		if ( chop( y ) == 0 && TY_W( x ) != 0 && TY_Q( x ) != 0 )
+		{
+			var_d1 = TY_X( x );
+			var_d2 = x;
+			
+			// solution of linear system for given d1, d2 to obtain a,b,ci,di
+			TY_SolveLinearEquations( var_d1, var_d2, &var_a, &var_b, &var_c1, &var_c2 );
+			
+			// select physical solutions, for details check paper: "Cluster formation in
+			// two-Yukawa fluids", J. Chem. Phys. 122 (2005)
+			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 )
+			{
+				sol_a[j]  = var_a;
+				sol_b[j]  = var_b;
+				sol_c1[j] = var_c1;
+				sol_c2[j] = var_c2;
+				sol_d1[j] = var_d1;
+				sol_d2[j] = var_d2;
+/*				
+				if ( debug )
+				{
+					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] ) );
+					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] ) );
+					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] ) );
+					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] ) );
+					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] ) );
+					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] ) );
+					
+					sprintf(buf, "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 );
+					XOPNotice(buf);
+				}
+ */
+				j++;
+			}
+		}
+	}
+	// number  remaining roots 
+	n_roots = j;
+	
+	
+	// if there is still more than one root left, than choose the one with the minimum
+	// average value inside the hardcore
+	if ( n_roots > 1 )
+	{
+		// the number of q values should be a power of 2
+		// in order to speed up the FFT
+		n = 1 << 14;
+		
+		// the maximum q value should be large enough 
+		// to enable a reasoble approximation of g(r)
+		qmax = 16 * 10 * 2 * Pi;
+		dq = qmax / ( n - 1 );
+		
+		// step size for g(r) = dr
+		
+		// allocate memory for pair correlation function g(r)
+		// and structure factor S(q)
+		sq = malloc( sizeof( double ) * n );
+		gr = malloc( sizeof( double ) * n );
+		
+		// loop over all remaining roots
+		min = 1e99;
+		selected_root = 10;	
+		sum = 0;
+		for ( j = 0; j < n_roots; j++) 
+		{
+			// calculate structure factor at different q values
+			for ( i = 0; i < n; i++) 
+			{
+				q = dq * i;
+				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] );
+/*				
+				if(i<20 && debug) {
+					sprintf(buf, "after SqTwoYukawa: s(q) = %g\r",sq[i] );
+					XOPNotice(buf);	
+				}
+ */
+			}
+			
+			// calculate pair correlation function for given
+			// structure factor, g(r) is computed at values
+			// r(i) = i * dr
+			PairCorrelation( phi, dq, sq, &dr, gr, n );
+			// determine sum inside the hardcore 
+			// 0 =< r < 1 of the pair-correlation function
+			sum = 0;
+			for (i = 0; i < floor( 1. / dr ); i++ )  {
+				sum += fabs( gr[i] );
+/*				
+				if(i<20 && debug) {
+					sprintf(buf, "g(r) in core = %g\r",fabs(gr[i]));
+					XOPNotice(buf);
+				}
+*/				
+			}
+
+			if ( sum < min )
+			{
+				min = sum;
+				selected_root = j;
+			}
+		}	
+		free( gr );
+		free( sq );
+		
+		// physical solution was found
+		*a  = sol_a [ selected_root ];//sol_a [ selected_root ];
+		*b  = sol_b [ selected_root ];
+		*c1 = sol_c1[ selected_root ];
+		*c2 = sol_c2[ selected_root ];
+		*d1 = sol_d1[ selected_root ];
+		*d2 = sol_d2[ selected_root ];
+		
+		return 1;
+	}
+	else if ( n_roots == 1 ) 
+	{
+		*a  = sol_a [0];
+		*b  = sol_b [0];
+		*c1 = sol_c1[0];
+		*c2 = sol_c2[0];
+		*d1 = sol_d1[0];
+		*d2 = sol_d2[0];
+		
+		return 1;
+	}
+
+	// no solution was found
+	return 0;
+}
Index: sansmodels/src/libigor/2Y_TwoYukawa.h
===================================================================
--- sansmodels/src/libigor/2Y_TwoYukawa.h	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
+++ sansmodels/src/libigor/2Y_TwoYukawa.h	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
@@ -0,0 +1,18 @@
+/*
+ *  TwoYukawa.h
+ *  TwoYukawa
+ *
+ *  Created by Marcus Hennig on 5/8/10.
+ *  Copyright 2010 __MyCompanyName__. All rights reserved.
+ *
+ */
+int TY_SolveEquations( double Z1, double Z2, double K1, double K2, double phi, 
+				       double a[], double b[], double c1[], double c2[], double d1[], double d2[], 
+				       int debug );
+
+int TY_CheckSolution( double Z1, double Z2, double K1, double K2, double phi, 
+				      double a, double b, double c1, double c2, double d1, double d2 );
+
+double SqTwoYukawa( double q, 
+			        double Z1, double Z2, double K1, double K2, double phi,
+			        double a, double b, double c1, double c2, double d1, double d2 );
Index: sansmodels/src/libigor/2Y_cpoly.c
===================================================================
--- sansmodels/src/libigor/2Y_cpoly.c	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
+++ sansmodels/src/libigor/2Y_cpoly.c	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
@@ -0,0 +1,660 @@
+/*
+ *******************************************************************************
+ *
+ *
+ *                       Copyright (c) 2002
+ *                       Henrik Vestermark
+ *                       Denmark
+ *
+ *                       All Rights Reserved
+ *
+ *   This source file is subject to the terms and conditions of the
+ *   Henrik Vestermark Software License Agreement which restricts the manner
+ *   in which it may be used.
+ *
+ *******************************************************************************
+ *
+ *
+ * Module name     :   cpoly.cpp
+ * Module ID Nbr   :   
+ * Description     :   cpoly.cpp -- Jenkins-Traub real polynomial root finder.
+ *                     Translation of TOMS493 from FORTRAN to C. This
+ *                     implementation of Jenkins-Traub partially adapts
+ *                     the original code to a C environment by restruction
+ *                     many of the 'goto' controls to better fit a block
+ *                     structured form. It also eliminates the global memory
+ *                     allocation in favor of local, dynamic memory management.
+ *
+ *                     The calling conventions are slightly modified to return
+ *                     the number of roots found as the function value.
+ *
+ *                     INPUT:
+ *                     opr - double precision vector of real coefficients in order of
+ *                          decreasing powers.
+ *                     opi - double precision vector of imaginary coefficients in order of
+ *                          decreasing powers.
+ *                     degree - integer degree of polynomial
+ *
+ *                     OUTPUT:
+ *                     zeror,zeroi - output double precision vectors of the
+ *                          real and imaginary parts of the zeros.
+ *                            to be consistent with rpoly.cpp the zeros is inthe index
+ *                            [0..max_degree-1]
+ *
+ *                     RETURN:
+ *                     returnval:   -1 if leading coefficient is zero, otherwise
+ *                          number of roots found. 
+ * --------------------------------------------------------------------------
+ * Change Record   :   
+ *
+ * Version	Author/Date		Description of changes
+ * -------  -----------		----------------------
+ * 01.01		HVE/021101		Initial release
+ *
+ * -- SRK changed name of *pi to *piw (and all occurances in this file) to avoid the conflict with pi in Apple's fp.h
+ *
+ * End of Change Record
+ * --------------------------------------------------------------------------
+ */
+
+
+/* define version string */
+
+#include <stdlib.h>
+#include <math.h>
+#include <float.h>
+
+static double sr, si, tr, ti, pvr, pvi, are, mre, eta, infin;
+static int nn;
+static double *pr, *piw, *hr, *hi, *qpr, *qpi, *qhr, *qhi, *shr, *shi; 
+
+static void noshft( const int l1 );
+static void fxshft( const int l2, double *zr, double *zi, int *conv );
+static void vrshft( const int l3, double *zr, double *zi, int *conv );
+static void calct( int *bol );
+static void nexth( const int bol );
+static void polyev( const int nn, const double sr, const double si, const double pr[], const double piw[], double qr[], double qi[], double *pvr, double *pvi );
+static double errev( const int nn, const double qr[], const double qi[], const double ms, const double mp, const double are, const double mre );
+static void cauchy( const int nn, double pt[], double q[], double *fn_val );
+static double scale( const int nn, const double pt[], const double eta, const double infin, const double smalno, const double base );
+static void cdivid( const double ar, const double ai, const double br, const double bi, double *cr, double *ci );
+static double cmod( const double r, const double i );
+static void mcon( double *eta, double *infiny, double *smalno, double *base );
+
+int cpoly( const double *opr, const double *opi, int degree, double *zeror, double *zeroi ) 
+{
+	int cnt1, cnt2, idnn2, i, conv;
+	double xx, yy, cosr, sinr, smalno, base, xxx, zr, zi, bnd;
+	
+	mcon( &eta, &infin, &smalno, &base );
+	are = eta;
+	mre = 2.0 * sqrt( 2.0 ) * eta;
+	xx = 0.70710678;
+	yy = -xx;
+	cosr = -0.060756474;
+	sinr = -0.99756405;
+	nn = degree;  
+	
+	// Algorithm fails if the leading coefficient is zero
+	if( opr[ 0 ] == 0 && opi[ 0 ] == 0 )
+		return -1;
+	
+	// Allocate arrays
+	pr = malloc( sizeof( double ) * ( degree + 1 ) );
+	piw = malloc( sizeof( double ) * ( degree + 1 ) );
+	hr =malloc( sizeof( double ) * ( degree + 1 ) );
+	hi = malloc( sizeof( double ) * ( degree + 1 ) );
+	qpr= malloc( sizeof( double ) * ( degree + 1 ) );
+	qpi= malloc( sizeof( double ) * ( degree + 1 ) );
+	qhr= malloc( sizeof( double ) * ( degree + 1 ) );
+	qhi= malloc( sizeof( double ) * ( degree + 1 ) );
+	shr= malloc( sizeof( double ) * ( degree + 1 ) );
+	shi= malloc( sizeof( double ) * ( degree + 1 ) );
+	
+	// Remove the zeros at the origin if any
+	while( opr[ nn ] == 0 && opi[ nn ] == 0 )
+	{
+		idnn2 = degree - nn;
+		zeror[ idnn2 ] = 0;
+		zeroi[ idnn2 ] = 0;
+		nn--;
+	}
+	
+	// Make a copy of the coefficients
+	for( i = 0; i <= nn; i++ )
+	{
+		pr[ i ] = opr[ i ];
+		piw[ i ] = opi[ i ];
+		shr[ i ] = cmod( pr[ i ], piw[ i ] );
+	}
+	
+	// Scale the polynomial
+	bnd = scale( nn, shr, eta, infin, smalno, base );
+	if( bnd != 1 )
+		for( i = 0; i <= nn; i++ )
+		{
+			pr[ i ] *= bnd;
+			piw[ i ] *= bnd;
+		}
+	
+search: 
+	if( nn <= 1 )
+	{
+		cdivid( -pr[ 1 ], -piw[ 1 ], pr[ 0 ], piw[ 0 ], &zeror[ degree-1 ], &zeroi[ degree-1 ] );
+		goto finish;
+	}
+	
+	// Calculate bnd, alower bound on the modulus of the zeros
+	for( i = 0; i<= nn; i++ )
+		shr[ i ] = cmod( pr[ i ], piw[ i ] );
+	
+	cauchy( nn, shr, shi, &bnd );
+	
+	// Outer loop to control 2 Major passes with different sequences of shifts
+	for( cnt1 = 1; cnt1 <= 2; cnt1++ )
+	{
+		// First stage  calculation , no shift
+		noshft( 5 );
+		
+		// Inner loop to select a shift
+		for( cnt2 = 1; cnt2 <= 9; cnt2++ )
+		{
+			// Shift is chosen with modulus bnd and amplitude rotated by 94 degree from the previous shif
+			xxx = cosr * xx - sinr * yy;
+			yy = sinr * xx + cosr * yy;
+			xx = xxx;
+			sr = bnd * xx;
+			si = bnd * yy;
+			
+			// Second stage calculation, fixed shift
+			fxshft( 10 * cnt2, &zr, &zi, &conv );
+			if( conv )
+            {
+				// The second stage jumps directly to the third stage ieration
+				// If successful the zero is stored and the polynomial deflated
+				idnn2 = degree - nn;
+				zeror[ idnn2 ] = zr;
+				zeroi[ idnn2 ] = zi;
+				nn--;
+				for( i = 0; i <= nn; i++ )
+				{
+					pr[ i ] = qpr[ i ];
+					piw[ i ] = qpi[ i ];
+				}
+				goto search;
+            }
+			// If the iteration is unsuccessful another shift is chosen
+		}
+		// if 9 shifts fail, the outer loop is repeated with another sequence of shifts
+	}
+	
+	// The zerofinder has failed on two major passes
+	// return empty handed with the number of roots found (less than the original degree)
+	degree -= nn;
+	
+finish:
+	// Deallocate arrays
+	free( pr );
+	free( piw );
+	free( hr );
+	free( hi );
+	free( qpr );
+	free( qpi );
+	free( qhr );
+	free( qhi );
+	free( shr );
+	free( shi );
+	
+	return degree;       
+}
+
+
+// COMPUTES  THE DERIVATIVE  POLYNOMIAL AS THE INITIAL H
+// POLYNOMIAL AND COMPUTES L1 NO-SHIFT H POLYNOMIALS.
+//
+static void noshft( const int l1 )
+{
+	int i, j, jj, n, nm1;
+	double xni, t1, t2;
+	
+	n = nn;
+	nm1 = n - 1;
+	for( i = 0; i < n; i++ )
+	{
+		xni = nn - i;
+		hr[ i ] = xni * pr[ i ] / n;
+		hi[ i ] = xni * piw[ i ] / n;
+	}
+	for( jj = 1; jj <= l1; jj++ )
+	{
+		if( cmod( hr[ n - 1 ], hi[ n - 1 ] ) > eta * 10 * cmod( pr[ n - 1 ], piw[ n - 1 ] ) )
+		{
+			cdivid( -pr[ nn ], -piw[ nn ], hr[ n - 1 ], hi[ n - 1 ], &tr, &ti );
+			for( i = 0; i < nm1; i++ )
+            {
+				j = nn - i - 1;
+				t1 = hr[ j - 1 ];
+				t2 = hi[ j - 1 ];
+				hr[ j ] = tr * t1 - ti * t2 + pr[ j ];
+				hi[ j ] = tr * t2 + ti * t1 + piw[ j ];
+            }
+			hr[ 0 ] = pr[ 0 ];
+			hi[ 0 ] = piw[ 0 ];
+		}
+		else
+		{
+			// If the constant term is essentially zero, shift H coefficients
+			for( i = 0; i < nm1; i++ )
+            {
+				j = nn - i - 1;
+				hr[ j ] = hr[ j - 1 ];
+				hi[ j ] = hi[ j - 1 ];
+            }
+			hr[ 0 ] = 0;
+			hi[ 0 ] = 0;
+		}
+	}
+}
+
+// COMPUTES L2 FIXED-SHIFT H POLYNOMIALS AND TESTS FOR CONVERGENCE.
+// INITIATES A VARIABLE-SHIFT ITERATION AND RETURNS WITH THE
+// APPROXIMATE ZERO IF SUCCESSFUL.
+// L2 - LIMIT OF FIXED SHIFT STEPS
+// ZR,ZI - APPROXIMATE ZERO IF CONV IS .TRUE.
+// CONV  - LOGICAL INDICATING CONVERGENCE OF STAGE 3 ITERATION
+//
+static void fxshft( const int l2, double *zr, double *zi, int *conv )
+{
+	int i, j, n;
+	int test, pasd, bol;
+	double otr, oti, svsr, svsi;
+	
+	n = nn;
+	polyev( nn, sr, si, pr, piw, qpr, qpi, &pvr, &pvi );
+	test = 1;
+	pasd = 0;
+	
+	// Calculate first T = -P(S)/H(S)
+	calct( &bol );
+	
+	// Main loop for second stage
+	for( j = 1; j <= l2; j++ )
+	{
+		otr = tr;
+		oti = ti;
+		
+		// Compute the next H Polynomial and new t
+		nexth( bol );
+		calct( &bol );
+		*zr = sr + tr;
+		*zi = si + ti;
+		
+		// Test for convergence unless stage 3 has failed once or this
+		// is the last H Polynomial
+		if( !( bol || !test || j == 12 ) )
+			if( cmod( tr - otr, ti - oti ) < 0.5 * cmod( *zr, *zi ) )
+            {
+				if( pasd )
+				{
+					// The weak convergence test has been passwed twice, start the third stage
+					// Iteration, after saving the current H polynomial and shift
+					for( i = 0; i < n; i++ )
+					{
+						shr[ i ] = hr[ i ];
+						shi[ i ] = hi[ i ];
+					}
+					svsr = sr;
+					svsi = si;
+					vrshft( 10, zr, zi, conv );
+					if( *conv ) return;
+					
+					//The iteration failed to converge. Turn off testing and restore h,s,pv and T
+					test = 0;
+					for( i = 0; i < n; i++ )
+					{
+						hr[ i ] = shr[ i ];
+						hi[ i ] = shi[ i ];
+					}
+					sr = svsr;
+					si = svsi;
+					polyev( nn, sr, si, pr, piw, qpr, qpi, &pvr, &pvi );
+					calct( &bol );
+					continue;
+				}
+				pasd = 1;
+            }
+			else
+				pasd = 0;
+	}
+	
+	// Attempt an iteration with final H polynomial from second stage
+	vrshft( 10, zr, zi, conv );
+}
+
+// CARRIES OUT THE THIRD STAGE ITERATION.
+// L3 - LIMIT OF STEPS IN STAGE 3.
+// ZR,ZI   - ON ENTRY CONTAINS THE INITIAL ITERATE, IF THE
+//           ITERATION CONVERGES IT CONTAINS THE FINAL ITERATE ON EXIT.
+// CONV    -  .TRUE. IF ITERATION CONVERGES
+//
+static void vrshft( const int l3, double *zr, double *zi, int *conv )
+{
+	int b, bol;
+	int i, j;
+	double mp, ms, omp, relstp, r1, r2, tp;
+	
+	*conv = 0;
+	b = 0;
+	sr = *zr;
+	si = *zi;
+	
+	// Main loop for stage three
+	for( i = 1; i <= l3; i++ )
+	{
+		// Evaluate P at S and test for convergence
+		polyev( nn, sr, si, pr, piw, qpr, qpi, &pvr, &pvi );
+		mp = cmod( pvr, pvi );
+		ms = cmod( sr, si );
+		if( mp <= 20 * errev( nn, qpr, qpi, ms, mp, are, mre ) )
+		{
+			// Polynomial value is smaller in value than a bound onthe error
+			// in evaluationg P, terminate the ietartion
+			*conv = 1;
+			*zr = sr;
+			*zi = si;
+			return;
+		}
+		if( i != 1 )
+		{
+			if( !( b || mp < omp || relstp >= 0.05 ) )
+            {
+				// Iteration has stalled. Probably a cluster of zeros. Do 5 fixed 
+				// shift steps into the cluster to force one zero to dominate
+				tp = relstp;
+				b = 1;
+				if( relstp < eta ) tp = eta;
+				r1 = sqrt( tp );
+				r2 = sr * ( 1 + r1 ) - si * r1;
+				si = sr * r1 + si * ( 1 + r1 );
+				sr = r2;
+				polyev( nn, sr, si, pr, piw, qpr, qpi, &pvr, &pvi );
+				for( j = 1; j <= 5; j++ )
+				{
+					calct( &bol );
+					nexth( bol );
+				}
+				omp = infin;
+				goto _20;
+            }
+			
+			// Exit if polynomial value increase significantly
+			if( mp *0.1 > omp ) return;
+		}
+		
+		omp = mp;
+		
+		// Calculate next iterate
+	_20:  calct( &bol );
+		nexth( bol );
+		calct( &bol );
+		if( !bol )
+		{
+			relstp = cmod( tr, ti ) / cmod( sr, si );
+			sr += tr;
+			si += ti;
+		}
+	}
+}
+
+// COMPUTES  T = -P(S)/H(S).
+// BOOL   - LOGICAL, SET TRUE IF H(S) IS ESSENTIALLY ZERO.
+static void calct( int *bol )
+{
+	int n;
+	double hvr, hvi;
+	
+	n = nn;
+	
+	// evaluate h(s)
+	polyev( n - 1, sr, si, hr, hi, qhr, qhi, &hvr, &hvi );
+	*bol = cmod( hvr, hvi ) <= are * 10 * cmod( hr[ n - 1 ], hi[ n - 1 ] ) ? 1 : 0;
+	if( !*bol )
+	{
+		cdivid( -pvr, -pvi, hvr, hvi, &tr, &ti );
+		return;
+	}
+	
+	tr = 0;
+	ti = 0;
+}
+
+// CALCULATES THE NEXT SHIFTED H POLYNOMIAL.
+// BOOL   -  LOGICAL, IF .TRUE. H(S) IS ESSENTIALLY ZERO
+//
+static void nexth( const int bol )
+{
+	int j, n;
+	double t1, t2;
+	
+	n = nn;
+	if( !bol )
+	{
+		for( j = 1; j < n; j++ )
+		{
+			t1 = qhr[ j - 1 ];
+			t2 = qhi[ j - 1 ];
+			hr[ j ] = tr * t1 - ti * t2 + qpr[ j ];
+			hi[ j ] = tr * t2 + ti * t1 + qpi[ j ];
+		}
+		hr[ 0 ] = qpr[ 0 ];
+		hi[ 0 ] = qpi[ 0 ];
+		return;
+	}
+	
+	// If h[s] is zero replace H with qh
+	for( j = 1; j < n; j++ )
+	{
+		hr[ j ] = qhr[ j - 1 ];
+		hi[ j ] = qhi[ j - 1 ];
+	}
+	hr[ 0 ] = 0;
+	hi[ 0 ] = 0;
+}
+
+// EVALUATES A POLYNOMIAL  P  AT  S  BY THE HORNER RECURRENCE
+// PLACING THE PARTIAL SUMS IN Q AND THE COMPUTED VALUE IN PV.
+//  
+static void polyev( const int nn, const double sr, const double si, const double pr[], const double piw[], double qr[], double qi[], double *pvr, double *pvi )  
+{
+	int i;
+	double t;
+	
+	qr[ 0 ] = pr[ 0 ];
+	qi[ 0 ] = piw[ 0 ];
+	*pvr = qr[ 0 ];
+	*pvi = qi[ 0 ];
+	
+	for( i = 1; i <= nn; i++ )
+	{
+		t = ( *pvr ) * sr - ( *pvi ) * si + pr[ i ];
+		*pvi = ( *pvr ) * si + ( *pvi ) * sr + piw[ i ];
+		*pvr = t;
+		qr[ i ] = *pvr;
+		qi[ i ] = *pvi;
+	}
+}
+
+// BOUNDS THE ERROR IN EVALUATING THE POLYNOMIAL BY THE HORNER RECURRENCE.
+// QR,QI - THE PARTIAL SUMS
+// MS    -MODULUS OF THE POINT
+// MP    -MODULUS OF POLYNOMIAL VALUE
+// ARE, MRE -ERROR BOUNDS ON COMPLEX ADDITION AND MULTIPLICATION
+//
+static double errev( const int nn, const double qr[], const double qi[], const double ms, const double mp, const double are, const double mre )
+{
+	int i;
+	double e;
+	
+	e = cmod( qr[ 0 ], qi[ 0 ] ) * mre / ( are + mre );
+	for( i = 0; i <= nn; i++ )
+		e = e * ms + cmod( qr[ i ], qi[ i ] );
+	
+	return e * ( are + mre ) - mp * mre;
+}
+
+// CAUCHY COMPUTES A LOWER BOUND ON THE MODULI OF THE ZEROS OF A
+// POLYNOMIAL - PT IS THE MODULUS OF THE COEFFICIENTS.
+//
+static void cauchy( const int nn, double pt[], double q[], double *fn_val )
+{
+	int i, n;
+	double x, xm, f, dx, df;
+	
+	pt[ nn ] = -pt[ nn ];
+	
+	// Compute upper estimate bound
+	n = nn;
+	x = exp( log( -pt[ nn ] ) - log( pt[ 0 ] ) ) / n;
+	if( pt[ n - 1 ] != 0 )
+	{
+		// Newton step at the origin is better, use it
+		xm = -pt[ nn ] / pt[ n - 1 ];
+		if( xm < x ) x = xm;
+	}
+	
+	// Chop the interval (0,x) until f < 0
+	while(1)
+	{
+		xm = x * 0.1;
+		f = pt[ 0 ];
+		for( i = 1; i <= nn; i++ )
+			f = f * xm + pt[ i ];
+		if( f <= 0 )
+			break;
+		x = xm;
+	}
+	dx = x;
+	
+	// Do Newton iteration until x converges to two decimal places
+	while( fabs( dx / x ) > 0.005 )
+	{
+		q[ 0 ] = pt[ 0 ];
+		for( i = 1; i <= nn; i++ )
+			q[ i ] = q[ i - 1 ] * x + pt[ i ];
+		f = q[ nn ];
+		df = q[ 0 ];
+		for( i = 1; i < n; i++ )
+			df = df * x + q[ i ];
+		dx = f / df;
+		x -= dx;
+	}
+	
+	*fn_val = x;
+}
+
+// RETURNS A SCALE FACTOR TO MULTIPLY THE COEFFICIENTS OF THE POLYNOMIAL.
+// THE SCALING IS DONE TO AVOID OVERFLOW AND TO AVOID UNDETECTED UNDERFLOW
+// INTERFERING WITH THE CONVERGENCE CRITERION.  THE FACTOR IS A POWER OF THE
+// BASE.
+// PT - MODULUS OF COEFFICIENTS OF P
+// ETA, INFIN, SMALNO, BASE - CONSTANTS DESCRIBING THE FLOATING POINT ARITHMETIC.
+//
+static double scale( const int nn, const double pt[], const double eta, const double infin, const double smalno, const double base )
+{
+	int i, l;
+	double hi, lo, max, min, x, sc;
+	double fn_val;
+	
+	// Find largest and smallest moduli of coefficients
+	hi = sqrt( infin );
+	lo = smalno / eta;
+	max = 0;
+	min = infin;
+	
+	for( i = 0; i <= nn; i++ )
+	{
+		x = pt[ i ];
+		if( x > max ) max = x;
+		if( x != 0 && x < min ) min = x;
+	}
+	
+	// Scale only if there are very large or very small components
+	fn_val = 1;
+	if( min >= lo && max <= hi ) return fn_val;
+	x = lo / min;
+	if( x <= 1 )
+		sc = 1 / ( sqrt( max )* sqrt( min ) );
+	else
+	{
+		sc = x;
+		if( infin / sc > max ) sc = 1;
+	}
+	l = (int)( log( sc ) / log(base ) + 0.5 );
+	fn_val = pow( base , l );
+	return fn_val;
+}
+
+// COMPLEX DIVISION C = A/B, AVOIDING OVERFLOW.
+//
+static void cdivid( const double ar, const double ai, const double br, const double bi, double *cr, double *ci )
+{
+	double r, d, t, infin;
+	
+	if( br == 0 && bi == 0 )
+	{
+		// Division by zero, c = infinity
+		mcon( &t, &infin, &t, &t );
+		*cr = infin;
+		*ci = infin;
+		return;
+	}
+	
+	if( fabs( br ) < fabs( bi ) )
+	{
+		r = br/ bi;
+		d = bi + r * br;
+		*cr = ( ar * r + ai ) / d;
+		*ci = ( ai * r - ar ) / d;
+		return;
+	}
+	
+	r = bi / br;
+	d = br + r * bi;
+	*cr = ( ar + ai * r ) / d;
+	*ci = ( ai - ar * r ) / d;
+}
+
+// MODULUS OF A COMPLEX NUMBER AVOIDING OVERFLOW.
+//
+static double cmod( const double r, const double i )
+{
+	double ar, ai;
+	
+	ar = fabs( r );
+	ai = fabs( i );
+	if( ar < ai )
+		return ai * sqrt( 1.0 + pow( ( ar / ai ), 2.0 ) );
+	
+	if( ar > ai )
+		return ar * sqrt( 1.0 + pow( ( ai / ar ), 2.0 ) );
+	
+	return ar * sqrt( 2.0 );
+}
+
+// MCON PROVIDES MACHINE CONSTANTS USED IN VARIOUS PARTS OF THE PROGRAM.
+// THE USER MAY EITHER SET THEM DIRECTLY OR USE THE STATEMENTS BELOW TO
+// COMPUTE THEM. THE MEANING OF THE FOUR CONSTANTS ARE -
+// ETA       THE MAXIMUM RELATIVE REPRESENTATION ERROR WHICH CAN BE DESCRIBED
+//           AS THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
+//           1.0_dp + ETA &gt; 1.0.
+// INFINY    THE LARGEST FLOATING-POINT NUMBER
+// SMALNO    THE SMALLEST POSITIVE FLOATING-POINT NUMBER
+// BASE      THE BASE OF THE FLOATING-POINT NUMBER SYSTEM USED
+//
+static void mcon( double *eta, double *infiny, double *smalno, double *base )
+{
+	*base = 10;//DBL_RADIX;
+	*eta = DBL_EPSILON;
+	*infiny = DBL_MAX;
+	*smalno = DBL_MIN;
+}
Index: sansmodels/src/libigor/2Y_cpoly.h
===================================================================
--- sansmodels/src/libigor/2Y_cpoly.h	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
+++ sansmodels/src/libigor/2Y_cpoly.h	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
@@ -0,0 +1,9 @@
+/*
+ *  cpoly.h
+ *  TwoYukawa
+ *
+ *  Created by Marcus Hennig on 5/8/10.
+ *  Copyright 2010 __MyCompanyName__. All rights reserved.
+ *
+ */
+int cpoly( const double *opr, const double *opi, int degree, double *zeror, double *zeroi ); 
Index: sansmodels/src/libigor/2Y_utility.c
===================================================================
--- sansmodels/src/libigor/2Y_utility.c	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
+++ sansmodels/src/libigor/2Y_utility.c	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
@@ -0,0 +1,17 @@
+/*
+ *  utility.c
+ *  twoyukawa
+ *
+ *  Created by Marcus Hennig on 5/13/10.
+ *  Copyright 2010 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "2Y_utility.h"
+double chop( double x )
+{
+	if ( fabs(x) < 1E-6 )
+		return 0;
+	else 
+		return x;
+}
Index: sansmodels/src/libigor/2Y_utility.h
===================================================================
--- sansmodels/src/libigor/2Y_utility.h	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
+++ sansmodels/src/libigor/2Y_utility.h	(revision a42ea1569853f2baa537de3a72362a8209d8b435)
@@ -0,0 +1,11 @@
+/*
+ *  utility.h
+ *  twoyukawa
+ *
+ *  Created by Marcus Hennig on 5/13/10.
+ *  Copyright 2010 __MyCompanyName__. All rights reserved.
+ *
+ */
+#include <math.h>
+
+double chop( double x );
