## finding those roots..

Discuss all kind of algorithms and data structures from their mathematical and programming sides.

Moderators: Darobat, RecursiveS, Dante Shamest, Bugdude, Wizard

### finding those roots..

Below I will try and explain the actual problem I am having.
I am trying to obtain a dxf file from the c++ code below to create horizontal profiles or slices through a 3d form.
With my little knowledge of c++ I have managed to do the following, which does work, however…
Code: Select all
`#include <fstream.h>#include <iostream.h>#include <math.h>#define   MaxNo 501void DXFSetUp(void);void DXFLines(void);void DXFText(void);void DXFFinishOff(void);int   i,j,NoSegments,Layer,Colour, hereornot[MaxNo];double x[MaxNo],y[MaxNo],z[MaxNo],         x3[MaxNo], y3[MaxNo], R[MaxNo], discriminant[MaxNo],        m1[MaxNo],m2[MaxNo],root[MaxNo],       xx1, yy1, xx2, yy2,        k, k1, k2,        alpha2, alpha3, lambda, mu, W,        v, e, u, a, PI,                root1, root2, root3, root4,       DXFx1,DXFy1,DXFz1,       DXFx2,DXFy2,DXFz2;        ofstream Chris("9complete_form.dxf");int main(void){DXFSetUp();xx1=0.05; yy1=0.05;xx2=0.03;yy2=0.9;alpha2=0.3; alpha3=1.5;mu=2.0;W=2.0;e=0.5;k=(sin(alpha2+alpha3))/(sin(alpha2)*sin(alpha2));k1=sin(alpha3)*sin(alpha3);k2=sin(alpha2)*sin(alpha2);NoSegments=MaxNo-1;for(v=0.1;v<=0.9;v+=0.01)//v=0.4;{PI=4.0*atan(1.0);for(i=0;i<=NoSegments;i+=1){root1=0.00001;  //these roots are for v=0.4root2=0.118946;root3=0.264714;root4=0.453145;u=root1+((root2-root1)/2)*(cos((i*PI)/NoSegments-PI)+1);            a=(k2*(1-v))/(k1*v+k2*(1-v));if ((a-e+u)!=0 ){x3[i]=(a*(e*xx2+u*(xx1-xx2)))/(e*(a-e+u));y3[i]=(a*(e*yy2+u*(yy1-yy2)))/(e*(a-e+u));hereornot[i]=1;}elsehereornot[i]=0;R[i]=(W*(e-u)/u)*pow(u/e,mu);discriminant[i]=(2*(xx1-x3[i])*(yy1-y3[i]))*(2*(xx1-x3[i])*(yy1-y3[i]))               -4*(R[i]*R[i]-(xx1-x3[i])*(xx1-x3[i]))*(R[i]*R[i]               -(yy1-y3[i])*(yy1-y3[i]));if(discriminant[i]>=0){m1[i]=((-2*(xx1-x3[i])*(yy1-y3[i]))+sqrt(discriminant[i]))/(2*(R[i]*R[i]-(xx1-x3[i])*(xx1-x3[i])));m2[i]=((-2*(xx1-x3[i])*(yy1-y3[i]))-sqrt(discriminant[i]))/(2*(R[i]*R[i]-(xx1-x3[i])*(xx1-x3[i])));hereornot[i]=1;}else hereornot[i]=0;if(a-e+u<0)root[i]=m1[i];elseroot[i]=m2[i];y[i]=(k*v*(y3[i]-root[i]*x3[i]))/(1+root[i]*root[i]);x[i]=-root[i]*y[i];z[i]=5*v;}Layer=0;Colour=10;for(i=1;i<=NoSegments;i+=1)if (hereornot[i-1]==1 && hereornot[i]==1){DXFx1=x[i-1];DXFy1=y[i-1];DXFz1=z[i-1];DXFx2=x[i];DXFy2=y[i];DXFz2=z[i];DXFLines();}for(j=0;j<=NoSegments;j+=1){u=root1+((root2-root1)/2)*(cos((j*PI)/NoSegments-PI)+1);if ((a-e+u)!=0 ){x3[j]=(a*(e*xx2+u*(xx1-xx2)))/(e*(a-e+u));y3[j]=(a*(e*yy2+u*(yy1-yy2)))/(e*(a-e+u));hereornot[j]=1;}elsehereornot[j]=0;R[j]=(W*(e-u)/u)*pow(u/e,mu);discriminant[j]=(2*(xx1-x3[j])*(yy1-y3[j]))*(2*(xx1-x3[j])*(yy1-y3[j]))               -4*(R[j]*R[j]-(xx1-x3[j])*(xx1-x3[j]))*(R[j]*R[j]               -(yy1-y3[j])*(yy1-y3[j]));if(discriminant[j]>=0){m1[j]=((-2*(xx1-x3[j])*(yy1-y3[j]))+sqrt(discriminant[j]))/(2*(R[j]*R[j]-(xx1-x3[j])*(xx1-x3[j])));m2[j]=((-2*(xx1-x3[j])*(yy1-y3[j]))-sqrt(discriminant[j]))/(2*(R[j]*R[j]-(xx1-x3[j])*(xx1-x3[j])));hereornot[j]=1;}else hereornot[j]=0;if(a-e+u<0)root[j]=m2[j];elseroot[j]=m1[j];y[j]=(k*v*(y3[j]-root[j]*x3[j]))/(1+root[j]*root[j]);x[j]=-root[j]*y[j];z[j]=5*v;}Layer=0;Colour=0;for(j=1;j<=NoSegments;j+=1)if (hereornot[j-1]==1 && hereornot[j]==1){DXFx1=x[j-1];DXFy1=y[j-1];DXFz1=z[j-1];DXFx2=x[j];DXFy2=y[j];DXFz2=z[j];DXFLines();}{if(root3>0&&root4>0)  {for(i=0;i<=NoSegments;i+=1){u=root3+((root4-root3)/2)*(cos((i*PI)/NoSegments-PI)+1);            a=(k2*(1-v))/(k1*v+k2*(1-v));if ((a-e+u)!=0 ){x3[i]=(a*(e*xx2+u*(xx1-xx2)))/(e*(a-e+u));y3[i]=(a*(e*yy2+u*(yy1-yy2)))/(e*(a-e+u));hereornot[i]=1;}elsehereornot[i]=0;R[i]=(W*(e-u)/u)*pow(u/e,mu);discriminant[i]=(2*(xx1-x3[i])*(yy1-y3[i]))*(2*(xx1-x3[i])*(yy1-y3[i]))               -4*(R[i]*R[i]-(xx1-x3[i])*(xx1-x3[i]))*(R[i]*R[i]               -(yy1-y3[i])*(yy1-y3[i]));if(discriminant[i]>=0){m1[i]=((-2*(xx1-x3[i])*(yy1-y3[i]))+sqrt(discriminant[i]))/(2*(R[i]*R[i]-(xx1-x3[i])*(xx1-x3[i])));m2[i]=((-2*(xx1-x3[i])*(yy1-y3[i]))-sqrt(discriminant[i]))/(2*(R[i]*R[i]-(xx1-x3[i])*(xx1-x3[i])));hereornot[i]=1;}else hereornot[i]=0;if(a-e+u<0)root[i]=m1[i];elseroot[i]=m2[i];y[i]=(k*v*(y3[i]-root[i]*x3[i]))/(1+root[i]*root[i]);x[i]=-root[i]*y[i];z[i]=5*v;}Layer=0;Colour=10;for(i=1;i<=NoSegments;i+=1)if (hereornot[i-1]==1 && hereornot[i]==1){DXFx1=x[i-1];DXFy1=y[i-1];DXFz1=z[i-1];DXFx2=x[i];DXFy2=y[i];DXFz2=z[i];DXFLines();}for(j=0;j<=NoSegments;j+=1){u=root3+((root4-root3)/2)*(cos((j*PI)/NoSegments-PI)+1);if ((a-e+u)!=0 ){x3[j]=(a*(e*xx2+u*(xx1-xx2)))/(e*(a-e+u));y3[j]=(a*(e*yy2+u*(yy1-yy2)))/(e*(a-e+u));hereornot[j]=1;}elsehereornot[j]=0;R[j]=(W*(e-u)/u)*pow(u/e,mu);discriminant[j]=(2*(xx1-x3[j])*(yy1-y3[j]))*(2*(xx1-x3[j])*(yy1-y3[j]))               -4*(R[j]*R[j]-(xx1-x3[j])*(xx1-x3[j]))*(R[j]*R[j]               -(yy1-y3[j])*(yy1-y3[j]));if(discriminant[j]>=0){m1[j]=((-2*(xx1-x3[j])*(yy1-y3[j]))+sqrt(discriminant[j]))/(2*(R[j]*R[j]-(xx1-x3[j])*(xx1-x3[j])));m2[j]=((-2*(xx1-x3[j])*(yy1-y3[j]))-sqrt(discriminant[j]))/(2*(R[j]*R[j]-(xx1-x3[j])*(xx1-x3[j])));hereornot[j]=1;}else hereornot[j]=0;if(a-e+u<0)root[j]=m2[j];elseroot[j]=m1[j];y[j]=(k*v*(y3[j]-root[j]*x3[j]))/(1+root[j]*root[j]);x[j]=-root[j]*y[j];z[j]=5*v;}Layer=0;Colour=0;for(j=1;j<=NoSegments;j+=1)if (hereornot[j-1]==1 && hereornot[j]==1){DXFx1=x[j-1];DXFy1=y[j-1];DXFz1=z[j-1];DXFx2=x[j];DXFy2=y[j];DXFz2=z[j];DXFLines();}}   }}DXFFinishOff();cout<<"DXF file written, end of program\n";               return 0;        }void DXFSetUp(void){Chris<<"0\n";Chris<<"SECTION\n";Chris<<"2\n";Chris<<"ENTITIES\n";}void DXFLines(void){Chris<<"0\nLINE\n8\n"<<Layer<<"\n";Chris<<"10\n"<<DXFx1<<"\n";Chris<<"20\n"<<DXFy1<<"\n";Chris<<"30\n"<<DXFz1<<"\n";Chris<<"11\n"<<DXFx2<<"\n";Chris<<"21\n"<<DXFy2<<"\n";Chris<<"31\n"<<DXFz2<<"\n";Chris<<"62\n"<<Colour<<"\n";}void DXFFinishOff(void){Chris<<"0\n";Chris<<"ENDSEC\n";Chris<<"0\n";Chris<<"EOF\n";Chris.close();}`

In the above code the main variable used is u to create the curve for each slice. The height of the slice is denoted by v, and then v is incremented to create the full form.
At the moment, for each value of v, I have to solve an equation to obtain the roots, root1, root2, root3 and root4, used in the program. I am using the following program, shown below, to find the roots, and then input these roots into the above program to obtain another profile. What I want to do is be able to solve the equation as part of the above program instead of having to manually input the roots.
Code: Select all
`#include <iostream> #include <stdlib.h> #include <math.h>float MyFunc (float); // The function bool isPositive(float); // Returns if MyFunc is +ve or -ve for passed value float v=0.4;float xx1=0.05,yy1=0.05,xx2=0.03,yy2=0.9,alpha2=0.3,alpha3=1.5,mu=2.0,W=2.0,e=0.5;double k=(sin(alpha2+alpha3))/(sin(alpha2)*sin(alpha3));double k1=sin(alpha3)*sin(alpha3);double k2=sin(alpha2)*sin(alpha2);double a=(k2*(1-v))/(k1*v+k2*(1-v));int main(){ bool intSign; float upperBound=2*e, lowerBound=0.00001, i, incr=0.0001; // Specify the bounds and the increment   i=lowerBound; //Start at lower end of interval do{ // do while below upper bound (i.e. for all roots in this range) intSign=isPositive(i); // Sign at lower bound of current interval     do{ // Search for change in sign         i+=incr;     }while(intSign==isPositive(i)&& i<upperBound);       if(intSign!=isPositive(i)){         cout << "Root in interval [" << i-incr << "," << i << "]" << endl;         cout << "root is "<<(i-incr+i)/2<<endl;        system("PAUSE");    } } while(i<upperBound);           return 0; } float MyFunc(float x) { // Function      return         ((2.0*(xx1-a*(e*xx2+x*(xx1-xx2))/(e*(a-e+x)))*                  (yy1-a*(e*yy2+x*(yy1-yy2))/(e*(a-e+x)))*                   2.0*(xx1-a*(e*xx2+x*(xx1-xx2))/(e*(a-e+x)))*(yy1-a*                   (e*yy2+x*(yy1-yy2))/(e*(a-e+x))))-                  4*((W*(e-x)/x*pow(x/e,mu))*(W*(e-x)/x*pow(x/e,mu))-                  (xx1-a*(e*xx2+x*(xx1-xx2))/(e*(a-e+x)))*(xx1-a*                  (e*xx2+x*(xx1-xx2))/(e*(a-e+x))))*                  ((W*(e-x)/x*pow(x/e,mu))*(W*(e-x)/x*pow(x/e,mu))-                  (yy1-a*(e*yy2+x*(yy1-yy2))/(e*(a-e+x)))*                  (yy1-a*(e*yy2+x*(yy1-yy2))/(e*(a-e+x)))));           } bool isPositive(float x) { // Test if +ve or -ve (false denotes -ve, true +ve)     if(MyFunc(x)>=0) return true;     else return false; }`

Two other points-
Root1 appears always to be zero for the constants used, although to avoid the problem of singularities I have taken it as 0.00001.
For the smaller values of v there are only two roots, 0 and one other. In the above program I have managed to say that if I manually input root3 and root4 as zero they are ignored.

Adapting the above program to achieve this is where I am getting very stuck…
Any ideas?
Many thanks
Dan P.

dont double post

Justin

Posts: 158
Joined: Tue Sep 30, 2003 10:07 am
Location: CA

Return to Algorithms & Data Structures

### Who is online

Users browsing this forum: No registered users and 1 guest