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..

Postby Dan P. » Tue Apr 13, 2004 7:30 am

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 501

void 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.4

root2=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;
}
else
hereornot[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];
else
root[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;
}
else
hereornot[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];
else
root[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;
}
else
hereornot[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];
else
root[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;
}
else
hereornot[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];
else
root[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.
 

Postby Justin » Sat Apr 17, 2004 9:25 pm

dont double post
User avatar
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 3 guests