draw5.cc
0000#include <math.h>
0001#include <stdio.h>
0002#include "draw5.h"
0003
0004
0005#define d2r 0.01745329252
0006#define r2d 57.295779513
0007
0008int outbyte(int num,char *buf);
0009int outstr(char *buf); /* for zero-terminated string */
0010int flush();
0011char msg[200];
0012
0013CThreeD thd;
0014int setting(double a[],int b[]);
0015int drawobj(int id,double tae[]);
0016
0017/*********************************************************************
0018
0019 DegreeToTAE
0020
0021 Put angles in degree into TransAngleEuler structure
0022
0023*********************************************************************/
0024void DegreeToTAE(double theta,double phi,double psi,TransAngleEuler *dest)
0025{
0026 dest->theta=d2r*theta;
0027 dest->phi=d2r*phi;
0028 dest->psi=d2r*psi;
0029}
0030
0031/*********************************************************************
0032
0033 DegreeToAE
0034
0035 Put angles in degree into AngleEuler structure
0036
0037*********************************************************************/
0038void DegreeToAE(double theta,double phi,double psi,AngleEuler *dest)
0039{
0040 dest->theta=d2r*theta;
0041 dest->phi=d2r*phi;
0042 dest->psi=d2r*psi;
0043}
0044
0045
0046/*********************************************************************
0047
0048 LinesToTAM
0049
0050 Construct TransAngleMatrix parameter from 2 Line3D parameters
0051
0052 a line segment on one of the coordinate axis
0053 b line segment on one of the coordinate axis
0054 w a and b are yzAxis, zxAxis or xyAxis
0055 dest TransAngleMatrix constructed
0056*********************************************************************/
0057void LinesToTAM(Line3d *a,Line3d *b,int w,TransAngleMatrix *dest)
0058{
0059 Point3d x,y,z;
0060 double d;
0061
0062 dest->x=a->sx;
0063 dest->y=a->sy;
0064 dest->z=a->sz;
0065
0066 switch (w){
0067 case yzAxis:
0068 NormLine3d(a,&y,&d);
0069 NormLine3d(b,&z,&d);
0070 x.x=y.y*z.z-y.z*z.y;
0071 x.y=y.z*z.x-y.x*z.z;
0072 x.z=y.x*z.y-y.y*z.x;
0073 break;
0074 case zxAxis:
0075 NormLine3d(a,&z,&d);
0076 NormLine3d(b,&x,&d);
0077 y.x=z.y*x.z-z.z*x.y;
0078 y.y=z.z*x.x-z.x*x.z;
0079 y.z=z.x*x.y-z.y*x.x;
0080 break;
0081 case xyAxis:
0082 NormLine3d(a,&x,&d);
0083 NormLine3d(b,&y,&d);
0084 z.x=x.y*y.z-x.z*y.y;
0085 z.y=x.z*y.x-x.x*y.z;
0086 z.z=x.x*y.y-x.y*y.x;
0087 break;
0088 }
0089 dest->c[0][0]=x.x;
0090 dest->c[1][0]=x.y;
0091 dest->c[2][0]=x.z;
0092 dest->c[0][1]=y.x;
0093 dest->c[1][1]=y.y;
0094 dest->c[2][1]=y.z;
0095 dest->c[0][2]=z.x;
0096 dest->c[1][2]=z.y;
0097 dest->c[2][2]=z.z;
0098}
0099
0100/*********************************************************************
0101
0102 LinesToTAE
0103
0104 Construct TransAngleEuler parameter from 2 Line3D parameters
0105
0106 a line segment on one of the coordinate axis
0107 b line segment on one of the coordinate axis
0108 w a and b are yzAxis, zxAxis or xyAxis
0109 dest TransAngleEuler constructed
0110*********************************************************************/
0111void LinesToTAE(Line3d *a,Line3d *b,int w,TransAngleEuler *dest)
0112{
0113 TransAngleMatrix tmp;
0114
0115 LinesToTAM(a,b,w,&tmp);
0116 dest->x=tmp.x;
0117 dest->y=tmp.y;
0118 dest->z=tmp.z;
0119
0120 if(tmp.c[0][2]==0.&&tmp.c[1][2]==0.){
0121 dest->theta=0.;
0122 dest->phi=0.;
0123 dest->psi=atan2(tmp.c[1][0],tmp.c[0][0]);
0124 }
0125 else{
0126 double x1,y1;
0127 dest->theta=acos(tmp.c[2][2]);
0128 dest->phi=atan2(tmp.c[1][2],tmp.c[0][2])+1.5707963268;
0129 x1=cos(dest->phi); // ascending node
0130 y1=sin(dest->phi);
0131
0132 // (y1*x.z-z1*x.y)*z.x+(z1*x.x-x1*x.z)*z.y+(x1*x.y-y1*x.x)*z.z
0133 // x1*x.x+y1*x.y+z1*x.z
0134 // dest->psi=atan2(y1*x.z*z.x-x1*x.z*z.y+(x1*x.y-y1*x.x)*z.z,x1*x.x+y1*x.y);
0135
0136 dest->psi=atan2(y1*tmp.c[2][0]*tmp.c[0][2]-x1*tmp.c[2][0]*tmp.c[2][1]
0137 +(x1*tmp.c[1][0]-y1*tmp.c[0][0])*tmp.c[2][2],
0138 x1*tmp.c[0][0]+y1*tmp.c[1][0]);
0139 }
0140}
0141
0142
0143/*********************************************************************
0144
0145 PointsToAM
0146
0147 Construct AngleMatrix parameter from 2 Point3D parameters
0148
0149 a point on one of the coordinate axis
0150 b point on one of the coordinate axis
0151 w a and b are yzAxis, zxAxis or xyAxis
0152 dest AngleMatrix constructed
0153*********************************************************************/
0154void PointsToAM(Point3d *a,Point3d *b,int w,AngleMatrix *dest)
0155{
0156 Point3d x,y,z;
0157 double d;
0158
0159
0160 switch (w){
0161 case yzAxis:
0162 NormPoint3d(a,&y,&d);
0163 NormPoint3d(b,&z,&d);
0164 x.x=y.y*z.z-y.z*z.y;
0165 x.y=y.z*z.x-y.x*z.z;
0166 x.z=y.x*z.y-y.y*z.x;
0167 break;
0168 case zxAxis:
0169 NormPoint3d(a,&z,&d);
0170 NormPoint3d(b,&x,&d);
0171 y.x=z.y*x.z-z.z*x.y;
0172 y.y=z.z*x.x-z.x*x.z;
0173 y.z=z.x*x.y-z.y*x.x;
0174 break;
0175 case xyAxis:
0176 NormPoint3d(a,&x,&d);
0177 NormPoint3d(b,&y,&d);
0178 z.x=x.y*y.z-x.z*y.y;
0179 z.y=x.z*y.x-x.x*y.z;
0180 z.z=x.x*y.y-x.y*y.x;
0181 break;
0182 }
0183 dest->c[0][0]=x.x;
0184 dest->c[1][0]=x.y;
0185 dest->c[2][0]=x.z;
0186 dest->c[0][1]=y.x;
0187 dest->c[1][1]=y.y;
0188 dest->c[2][1]=y.z;
0189 dest->c[0][2]=z.x;
0190 dest->c[1][2]=z.y;
0191 dest->c[2][2]=z.z;
0192}
0193
0194/*********************************************************************
0195
0196 PointsToAE
0197
0198 Construct AngleEuler parameter from 2 Point3D parameters
0199
0200 a point on one of the coordinate axis
0201 b point on one of the coordinate axis
0202 w a and b are yzAxis, zxAxis or xyAxis
0203 dest TransAngleEuler constructed
0204*********************************************************************/
0205void PointsToAE(Point3d *a,Point3d *b,int w,AngleEuler *dest)
0206{
0207 AngleMatrix tmp;
0208
0209 PointsToAM(a,b,w,&tmp);
0210
0211 if(tmp.c[0][2]==0.&&tmp.c[1][2]==0.){
0212 dest->theta=0.;
0213 dest->phi=0.;
0214 dest->psi=atan2(tmp.c[1][0],tmp.c[0][0]);
0215 }
0216 else{
0217 double x1,y1;
0218 dest->theta=acos(tmp.c[2][2]);
0219 dest->phi=atan2(tmp.c[1][2],tmp.c[0][2])+1.5707963268;
0220 x1=cos(dest->phi); // ascending node
0221 y1=sin(dest->phi);
0222
0223 // (y1*x.z-z1*x.y)*z.x+(z1*x.x-x1*x.z)*z.y+(x1*x.y-y1*x.x)*z.z
0224 // x1*x.x+y1*x.y+z1*x.z
0225 // dest->psi=atan2(y1*x.z*z.x-x1*x.z*z.y+(x1*x.y-y1*x.x)*z.z,x1*x.x+y1*x.y);
0226
0227 dest->psi=atan2(y1*tmp.c[2][0]*tmp.c[0][2]-x1*tmp.c[2][0]*tmp.c[2][1]
0228 +(x1*tmp.c[1][0]-y1*tmp.c[0][0])*tmp.c[2][2],
0229 x1*tmp.c[0][0]+y1*tmp.c[1][0]);
0230 }
0231}
0232
0233
0234/*********************************************************************
0235 NormLine3d
0236
0237 normalize a vector
0238
0239 source line segment to be normalized
0240
0241 dest normalized vector
0242 norm length of the vector
0243
0244*********************************************************************/
0245void NormLine3d(Line3d *source,Point3d *dest,double *norm)
0246{
0247 double x,y,z,d;
0248
0249 x=source->ex-source->sx;
0250 y=source->ey-source->sy;
0251 z=source->ez-source->sz;
0252
0253 d=sqrt(x*x+y*y+z*z);
0254 dest->x=x/d;
0255 dest->y=y/d;
0256 dest->z=z/d;
0257 *norm=d;
0258}
0259
0260/*********************************************************************
0261 NormPoint3d
0262
0263 normalize a vector
0264
0265 source vector to be normalized
0266
0267 dest normalized vector
0268 norm length of the vector
0269
0270*********************************************************************/
0271void NormPoint3d(Point3d *source,Point3d *dest,double *norm)
0272{
0273 double x,y,z,d;
0274
0275 x=source->x;
0276 y=source->y;
0277 z=source->z;
0278
0279 d=sqrt(x*x+y*y+z*z);
0280 dest->x=x/d;
0281 dest->y=y/d;
0282 dest->z=z/d;
0283 *norm=d;
0284}
0285
0286/*********************************************************************
0287 NormLine2d
0288
0289 normalize a vector
0290
0291 source line segment to be normalized
0292
0293 dest normalized vector
0294 norm length of the vector
0295
0296*********************************************************************/
0297void NormLine2d(Line2d *source,Point2d *dest,double *norm)
0298{
0299 double x,y,d;
0300
0301 x=source->ex-source->sx;
0302 y=source->ey-source->sy;
0303
0304 d=sqrt(x*x+y*y);
0305 dest->x=x/d;
0306 dest->y=y/d;
0307 *norm=d;
0308}
0309
0310/*********************************************************************
0311 NormPoint2d
0312
0313 normalize a vector
0314
0315 source vector to be normalized
0316
0317 dest normalized vector
0318 norm length of the vector
0319
0320*********************************************************************/
0321void NormPoint2d(Point2d *source,Point2d *dest,double *norm)
0322{
0323 double x,y,d;
0324
0325 x=source->x;
0326 y=source->y;
0327
0328 d=sqrt(x*x+y*y);
0329 dest->x=x/d;
0330 dest->y=y/d;
0331 *norm=d;
0332}
0333
0334
0335/*********************************************************************
0336 TAEtoTAM
0337
0338 convert TransAngleEuler to TransAngleMatrix
0339
0340 source TransAngleEuler
0341
0342 dest TransAngleMatrix
0343
0344*********************************************************************/
0345void TAEtoTAM(TransAngleEuler *source,TransAngleMatrix *dest)
0346{
0347 double sth,cth,sph,cph,sps,cps;
0348
0349 dest->x=source->x;
0350 dest->y=source->y;
0351 dest->z=source->z;
0352 sth=sin(source->theta);
0353 cth=cos(source->theta);
0354 sph=sin(source->phi);
0355 cph=cos(source->phi);
0356 sps=sin(source->psi);
0357 cps=cos(source->psi);
0358 dest->c[0][0]=cps*cph-sps*cth*sph; // x component of x'
0359 dest->c[0][1]=-sps*cph-cps*cth*sph;
0360 dest->c[0][2]=sth*sph;
0361 dest->c[1][0]=cps*sph+sps*cth*cph; // y component of x'
0362 dest->c[1][1]=-sps*sph+cps*cth*cph;
0363 dest->c[1][2]=-sth*cph;
0364 dest->c[2][0]=sps*sth; // z component of x'
0365 dest->c[2][1]=cps*sth;
0366 dest->c[2][2]=cth;
0367}
0368
0369
0370/*********************************************************************
0371 TAMtoTAE
0372
0373 convert TransAngleMatrix to TransAngleEuler
0374
0375 source TransAngleMatrix
0376
0377 dest TransAngleEuler
0378
0379*********************************************************************/
0380void TAMtoTAE(TransAngleMatrix *source,TransAngleEuler *dest)
0381{
0382 double a,ca,sa;
0383
0384 dest->x=source->x;
0385 dest->y=source->y;
0386 dest->z=source->z;
0387
0388 if(source->c[2][2]==1.){
0389 dest->theta=0.;
0390 dest->phi=0.;
0391 dest->psi=atan2(source->c[0][1],source->c[0][0]);
0392 return;
0393 }
0394 dest->theta=acos(source->c[2][2]);
0395 dest->phi=atan2(source->c[1][2],source->c[0][2])+1.5707963268;
0396 a=dest->phi;
0397 ca=cos(a);
0398 sa=sin(a);
0399 dest->psi=-atan2(source->c[0][1]*ca+source->c[1][1]*sa,
0400 source->c[0][0]*ca+source->c[1][0]*sa);
0401}
0402
0403
0404/*********************************************************************
0405 InverseTAM
0406
0407 Inverse body-frame order
0408
0409 source TransAngleMatrix
0410
0411 dest TransAngleMatrix
0412
0413*********************************************************************/
0414void InverseTAM(TransAngleMatrix *source,TransAngleMatrix *dest)
0415{
0416 Point3d a;
0417
0418 a.x=-source->x;
0419 a.y=-source->y;
0420 a.z=-source->z;
0421
0422 dest->x=source->c[0][0]*a.x+source->c[1][0]*a.y+source->c[2][0]*a.z;
0423 dest->y=source->c[0][1]*a.x+source->c[1][1]*a.y+source->c[2][1]*a.z;
0424 dest->z=source->c[0][2]*a.x+source->c[1][2]*a.y+source->c[2][2]*a.z;
0425
0426 dest->c[0][0]=source->c[0][0];
0427 dest->c[0][1]=source->c[1][0];
0428 dest->c[0][2]=source->c[2][0];
0429 dest->c[1][0]=source->c[0][1];
0430 dest->c[1][1]=source->c[1][1];
0431 dest->c[1][2]=source->c[2][1];
0432 dest->c[2][0]=source->c[0][2];
0433 dest->c[2][1]=source->c[1][2];
0434 dest->c[2][2]=source->c[2][2];
0435
0436}
0437
0438
0439
0440/*********************************************************************
0441 BodyToFrameTAM
0442
0443 convert Body coordinate to Frame coordinate
0444
0445 source Point3d
0446 m TransAngleMatrix
0447
0448 dest Point3d
0449
0450*********************************************************************/
0451void BodyToFrameTAM(Point3d *source,TransAngleMatrix *m,Point3d *dest)
0452{
0453 dest->x=m->x+m->c[0][0]*source->x+m->c[0][1]*source->y+m->c[0][2]*source->z;
0454 dest->y=m->y+m->c[1][0]*source->x+m->c[1][1]*source->y+m->c[1][2]*source->z;
0455 dest->z=m->z+m->c[2][0]*source->x+m->c[2][1]*source->y+m->c[2][2]*source->z;
0456}
0457
0458/*********************************************************************
0459 FrameToBodyTAM
0460
0461 convert Frame coordinate to Body coordinate
0462
0463 source Point3d
0464 m TransAngleMatrix
0465
0466 dest Point3d
0467
0468*********************************************************************/
0469void FrameToBodyTAM(Point3d *source,TransAngleMatrix *m,Point3d *dest)
0470{
0471 dest->x=m->c[0][0]*(source->x-m->x)
0472 +m->c[1][0]*(source->y-m->y)
0473 +m->c[2][0]*(source->z-m->z);
0474 dest->y=m->c[0][1]*(source->x-m->x)
0475 +m->c[1][1]*(source->y-m->y)
0476 +m->c[2][1]*(source->z-m->z);
0477 dest->z=m->c[0][2]*(source->x-m->x)
0478 +m->c[1][2]*(source->y-m->y)
0479 +m->c[2][2]*(source->z-m->z);
0480}
0481
0482
0483/*********************************************************************
0484 MoveTwiceTAE
0485
0486 obtain single conversion parameters from double consecutive
0487 coordinate conversion parameters
0488
0489 src1 TransAngleEuler cood1 - frame cood2 - body
0490 src2 TransAngleEuler cood2 - frame cood3 - body
0491
0492 dest TransAngleEuler cood1 - frame cood3 - body
0493
0494********************************************************************/
0495void MoveTwiceTAE(TransAngleEuler *src1,TransAngleEuler *src2,TransAngleEuler *dest)
0496{
0497 double x1,y1,z1,theta1,phi1,psi1;
0498 double x2,y2,z2,theta2,phi2,psi2;
0499 double x3,y3,z3,theta3,phi3,psi3;
0500
0501 double a[3][3],b[3][3],c[3],d[3],e[3];
0502 double cost,sint,cosf,sinf,cosp,sinp;
0503 bool locked;
0504
0505 x1=src1->x;
0506 y1=src1->y;
0507 z1=src1->z;
0508 theta1=src1->theta;
0509 phi1=src1->phi;
0510 psi1=src1->psi;
0511 x2=src2->x;
0512 y2=src2->y;
0513 z2=src2->z;
0514 theta2=src2->theta;
0515 phi2=src2->phi;
0516 psi2=src2->psi;
0517
0518 cost=cos(theta1);
0519 sint=sin(theta1);
0520 cosf=cos(phi1);
0521 sinf=sin(phi1);
0522 cosp=cos(psi1);
0523 sinp=sin(psi1);
0524 // if d[i] is a vector in cood2, c[j]=a[j][i]*d[i] is the vector in cood1
0525 a[0][0]=cosp*cosf-sinp*cost*sinf;
0526 a[0][1]=-sinp*cosf-cosp*cost*sinf;
0527 a[0][2]=sint*sinf;
0528
0529 a[1][0]=cosp*sinf+sinp*cost*cosf;
0530 a[1][1]=-sinp*sinf+cosp*cost*cosf;
0531 a[1][2]=-sint*cosf;
0532
0533 a[2][0]=sinp*sint;
0534 a[2][1]=cosp*sint;
0535 a[2][2]=cost;
0536
0537 cost=cos(theta2);
0538 sint=sin(theta2);
0539 cosf=cos(phi2);
0540 sinf=sin(phi2);
0541 cosp=cos(psi2);
0542 sinp=sin(psi2);
0543
0544 // if d[i] is a vector in cood3, c[j]=b[j][i]*d[i] is the vector in cood2
0545 b[0][0]=cosp*cosf-sinp*cost*sinf;
0546 b[0][1]=-sinp*cosf-cosp*cost*sinf;
0547 b[0][2]=sint*sinf;
0548
0549 b[1][0]=cosp*sinf+sinp*cost*cosf;
0550 b[1][1]=-sinp*sinf+cosp*cost*cosf;
0551 b[1][2]=-sint*cosf;
0552
0553 b[2][0]=sinp*sint;
0554 b[2][1]=cosp*sint;
0555 b[2][2]=cost;
0556
0557 x3=x1+a[0][0]*x2+a[0][1]*y2+a[0][2]*z2;
0558 y3=y1+a[1][0]*x2+a[1][1]*y2+a[1][2]*z2;
0559 z3=z1+a[2][0]*x2+a[2][1]*y2+a[2][2]*z2;
0560
0561 // if c is a unit z vector in coord 3, d=b*c is the vector in coord 2
0562 // and a*d is the vector in coord 1
0563 c[0]=0.;
0564 c[1]=0.;
0565 c[2]=1.;
0566 d[0]=b[0][0]*c[0]+b[0][1]*c[1]+b[0][2]*c[2];
0567 d[1]=b[1][0]*c[0]+b[1][1]*c[1]+b[1][2]*c[2];
0568 d[2]=b[2][0]*c[0]+b[2][1]*c[1]+b[2][2]*c[2];
0569 e[0]=a[0][0]*d[0]+a[0][1]*d[1]+a[0][2]*d[2];
0570 e[1]=a[1][0]*d[0]+a[1][1]*d[1]+a[1][2]*d[2];
0571 e[2]=a[2][0]*d[0]+a[2][1]*d[1]+a[2][2]*d[2];
0572
0573 // Euler angles of coord3 with reference to coord1 are to be derived
0574 // theta is the angle between z vector in coord3 and z vector in coord1
0575 // cosine of the theta is e[2] in the above expression
0576
0577 theta3=acos(e[2]);
0578 if(theta3<0.0005) {
0579 phi3=0.; // the case when theta = 0
0580 locked=true;
0581 }
0582 else {
0583 phi3=atan2(e[0],-e[1]);
0584 locked=false;
0585 }
0586
0587 // if c is a unit x vector in coord 3, d=b*c is the vector in coord 2
0588 // and a*d is the vector in coord 1
0589 c[0]=1.;
0590 c[1]=0.;
0591 c[2]=0.;
0592 d[0]=b[0][0]*c[0]+b[0][1]*c[1]+b[0][2]*c[2];
0593 d[1]=b[1][0]*c[0]+b[1][1]*c[1]+b[1][2]*c[2];
0594 d[2]=b[2][0]*c[0]+b[2][1]*c[1]+b[2][2]*c[2];
0595 e[0]=a[0][0]*d[0]+a[0][1]*d[1]+a[0][2]*d[2];
0596 e[1]=a[1][0]*d[0]+a[1][1]*d[1]+a[1][2]*d[2];
0597 e[2]=a[2][0]*d[0]+a[2][1]*d[1]+a[2][2]*d[2];
0598 // if c is a unit y vector in coord 3, d=b*c is the vector in coord 2
0599 // and a*d is the vector in coord 1
0600 c[0]=0.;
0601 c[1]=1.;
0602 c[2]=0.;
0603 d[0]=b[0][0]*c[0]+b[0][1]*c[1]+b[0][2]*c[2];
0604 d[1]=b[1][0]*c[0]+b[1][1]*c[1]+b[1][2]*c[2];
0605 d[2]=b[2][0]*c[0]+b[2][1]*c[1]+b[2][2]*c[2];
0606 c[0]=a[0][0]*d[0]+a[0][1]*d[1]+a[0][2]*d[2];
0607 c[1]=a[1][0]*d[0]+a[1][1]*d[1]+a[1][2]*d[2];
0608 c[2]=a[2][0]*d[0]+a[2][1]*d[1]+a[2][2]*d[2];
0609
0610 if(locked) psi3=atan2(e[1],e[0]);
0611 else psi3=atan2(e[2],c[2]);
0612
0613 dest->x=x3;
0614 dest->y=y3;
0615 dest->z=z3;
0616 dest->theta=theta3;
0617 dest->phi=phi3;
0618 dest->psi=psi3;
0619}
0620
0621
0622/*********************************************************************
0623 AEtoAM
0624
0625 convert CoordEuler to CoordMatrix
0626
0627 source AngleEuler
0628
0629 dest AngleMatrix
0630
0631*********************************************************************/
0632void AEtoAM(AngleEuler *source,AngleMatrix *dest)
0633{
0634 double sth,cth,sph,cph,sps,cps;
0635
0636 sth=sin(source->theta);
0637 cth=cos(source->theta);
0638 sph=sin(source->phi);
0639 cph=cos(source->phi);
0640 sps=sin(source->psi);
0641 cps=cos(source->psi);
0642 dest->c[0][0]=cps*cph-sps*cth*sph;
0643 dest->c[0][1]=-sps*cph-cps*cth*sph;
0644 dest->c[0][2]=sth*sph;
0645 dest->c[1][0]=cps*sph+sps*cth*cph;
0646 dest->c[1][1]=-sps*sph+cps*cth*cph;
0647 dest->c[1][2]=-sth*cph;
0648 dest->c[2][0]=sps*sth;
0649 dest->c[2][1]=cps*sth;
0650 dest->c[2][2]=cth;
0651}
0652
0653/*********************************************************************
0654 BodyToFrameAM
0655
0656 convert Body coordinate to Frame coordinate
0657
0658 source Point3d - Cartesian
0659 m AngleMatrix
0660
0661 dest Point3d - Cartesian
0662
0663*********************************************************************/
0664void BodyToFrameAM(Point3d *source,AngleMatrix *m,Point3d *dest)
0665{
0666
0667 dest->x=m->c[0][0]*source->x+m->c[0][1]*source->y+m->c[0][2]*source->z;
0668 dest->y=m->c[1][0]*source->x+m->c[1][1]*source->y+m->c[1][2]*source->z;
0669 dest->z=m->c[2][0]*source->x+m->c[2][1]*source->y+m->c[2][2]*source->z;
0670
0671}
0672
0673/*********************************************************************
0674 FrameToBodyAM
0675
0676 convert Frame coordinate to Body coordinate
0677
0678 source Point3d - Cartesian
0679 m AngleMatrix
0680
0681 dest Point3d - Cartesian
0682
0683*********************************************************************/
0684void FrameToBodyAM(Point3d *source,AngleMatrix *m,Point3d *dest)
0685{
0686 dest->x=m->c[0][0]*source->x+m->c[1][0]*source->y+m->c[2][0]*source->z;
0687 dest->y=m->c[0][1]*source->x+m->c[1][1]*source->y+m->c[2][1]*source->z;
0688 dest->z=m->c[0][2]*source->x+m->c[1][2]*source->y+m->c[2][2]*source->z;
0689}
0690
0691/*********************************************************************
0692 BodyToFrameAMP
0693
0694 convert Body coordinate to Frame coordinate - Polar
0695
0696 source Point2d - azimuth and elevation (rad)
0697 m AngleMatrix
0698
0699 dest Point2d - azimuth and elevation (rad)
0700
0701*********************************************************************/
0702void BodyToFrameAMP(Point2d *source,AngleMatrix *m,Point2d *dest)
0703{
0704 Point3d s,d; // unit vector
0705
0706 ToCartesian3d(source,&s);
0707
0708 d.x=m->c[0][0]*s.x+m->c[0][1]*s.y+m->c[0][2]*s.z;
0709 d.y=m->c[1][0]*s.x+m->c[1][1]*s.y+m->c[1][2]*s.z;
0710 d.z=m->c[2][0]*s.x+m->c[2][1]*s.y+m->c[2][2]*s.z;
0711
0712 ToPolar3d(&d,dest);
0713
0714}
0715
0716/*********************************************************************
0717 FrameToBodyAMP
0718
0719 convert Frame coordinate to Body coordinate - Polar
0720
0721 source Point2d - azimuth and elevation (rad)
0722 m AngleMatrix
0723
0724 dest Point2d - azimuth and elevation (rad)
0725
0726*********************************************************************/
0727void FrameToBodyAMP(Point2d *source,AngleMatrix *m,Point2d *dest)
0728{
0729 Point3d s,d;
0730
0731 ToCartesian3d(source,&s);
0732
0733 d.x=m->c[0][0]*s.x+m->c[1][0]*s.y+m->c[2][0]*s.z;
0734 d.y=m->c[0][1]*s.x+m->c[1][1]*s.y+m->c[2][1]*s.z;
0735 d.z=m->c[0][2]*s.x+m->c[1][2]*s.y+m->c[2][2]*s.z;
0736
0737 ToPolar3d(&d,dest);
0738
0739}
0740
0741
0742/*********************************************************************
0743 ToPolar3d
0744
0745 convert Catesian coordinate to Polar coordinate
0746
0747 source Point3dDouble - Cartesian with unit length
0748
0749 dest Point2dDouble - azimuth and elevation (rad)
0750
0751*********************************************************************/
0752void ToPolar3d(Point3d *source,Point2d *dest)
0753{
0754 dest->y=asin(source->z);
0755 if(source->z==1.||source->z==-1.) dest->x=0.;
0756 else dest->x=atan2(source->y,source->x);
0757}
0758
0759/*********************************************************************
0760 ToCartesian3d
0761
0762 convert Polar coordinate to Cartesian coordinate
0763
0764 source Point2dDouble - azimuth and elevation (rad)
0765
0766 dest Point3dDouble - Cartesian with unit length
0767
0768*********************************************************************/
0769void ToCartesian3d(Point2d *source,Point3d *dest)
0770{
0771 double cel;
0772 cel=cos(source->y); // cosine elevation
0773 dest->x=cel*cos(source->x);
0774 dest->y=cel*sin(source->x);
0775 dest->z=sin(source->y);
0776}
0777
0778/*********************************************************************
0779 VectorProduct
0780
0781 Vector product
0782
0783 u input
0784 v input
0785 w output uxv
0786
0787*********************************************************************/
0788void VectorProduct(Point3d *u,Point3d *v,Point3d *w)
0789{
0790 w->x=u->y*v->z-u->z*v->y;
0791 w->y=u->z*v->x-u->x*v->z;
0792 w->z=u->x*v->y-u->y*v->x;
0793}
0794
0795/*********************************************************************
0796 AcuteAngle
0797
0798 Test if the angle u-v-w is an acute angle
0799
0800 u point1
0801 v point2 - center point
0802 w point3
0803
0804 returned true if the angle is acute angle
0805 false if the angle is obtuse or right angle
0806
0807*********************************************************************/
0808bool AcuteAngle(Point3d *u,Point3d *v,Point3d *w)
0809{
0810 if((u->x-v->x)*(w->x-v->x)+(u->y-v->y)*(w->y-v->y)+
0811 (u->z-v->z)*(w->z-v->z)>0.) return true;
0812 else return false;
0813}
0814
0815/**************************************************************************
0816 MakeThePolar
0817
0818 'polar' is a coordinate system with its origin at the viewer's eye.
0819 Its z axis is directed toward line of sight and x-z plane is pararel
0820 to the world z axis. This coordinate is for generating polar coordinate
0821 where the north pole is directed to the line of sight and the azimuth
0822 is measured crockwise from the direction pointing upward.
0823
0824 sight starts at the eye and ends at the target
0825
0826 polar with origin at the eye and z-axis toward the line of sight
0827**************************************************************************/
0828void MakeThePolar(Line3d *sight,TransAngleMatrix *polar)
0829{
0830 Point3d x,y,z;
0831 double d;
0832
0833 polar->x=sight->sx;
0834 polar->y=sight->sy;
0835 polar->z=sight->sz;
0836
0837 NormLine3d(sight,&z,&d);
0838
0839 // z.x z.y z.z
0840 // 0 0 1
0841
0842 y.x=z.y;
0843 y.y=-z.x;
0844 y.z=0.;
0845 NormPoint3d(&y,&y,&d);
0846
0847 // z.y -z.x 0
0848 // z.x z.y z.z
0849
0850 x.x=-z.x*z.z;
0851 x.y=-z.y*z.z;
0852 x.z=z.y*z.y+z.x*z.x;
0853 NormPoint3d(&x,&x,&d);
0854
0855 polar->c[0][0]=x.x;
0856 polar->c[0][1]=y.x;
0857 polar->c[0][2]=z.x;
0858 polar->c[1][0]=x.y;
0859 polar->c[1][1]=y.y;
0860 polar->c[1][2]=z.y;
0861 polar->c[2][0]=x.z;
0862 polar->c[2][1]=y.z;
0863 polar->c[2][2]=z.z;
0864}
0865
0866/**************************************************************************
0867 WorldToPolarCnv
0868
0869 'polar' is a coordinate system with its origin at the viewer's eye.
0870 Its z axis is directed toward line of sight and x-z plane is pararel
0871 to the world z axis. This coordinate is for generating polar coordinate
0872 where the north pole is directed to the line of sight and the azimuth
0873 is measured crockwise from the direction pointing upward.
0874
0875 polar with origin at the eye and z-axis toward the line of sight
0876 s point in World coordinate
0877
0878 d polar coordinate az(rad), el(rad), and distance
0879 el measured fron boresight
0880
0881**************************************************************************/
0882void WorldToPolarCnv(Point3d *s,TransAngleMatrix *polar,Point3d *d)
0883{
0884 Point3d a;
0885 double b;
0886
0887 FrameToBodyTAM(s,polar,d);
0888 NormPoint3d(d,&a,&b);
0889 if(a.z==1.){
0890 d->x=0.;
0891 d->y=0.;
0892 d->z=b;
0893 }
0894 else {
0895 d->x=atan2(a.y,a.x);
0896 d->y=acos(a.z);
0897 d->z=b;
0898 }
0899}
0900
0901/**************************************************************************
0902 WorldToPerspCnv
0903
0904 'polar' is a coordinate system with its origin at the viewer's eye.
0905 Its z axis is directed toward line of sight and x-z plane is pararel
0906 to the world z axis. This coordinate is for generating polar coordinate
0907 where the north pole is directed to the line of sight and the azimuth
0908 is measured crockwise from the direction pointing upward.
0909
0910 polar with origin at the eye and z-axis toward the line of sight
0911 s point in World coordinate
0912
0913 d polar coordinate az(rad), el(rad), and distance
0914 el measured fron boresight
0915
0916**************************************************************************/
0917void WorldToPerspCnv(Point3d *s,TransAngleMatrix *polar,PointPerspectve *d)
0918{
0919 Point3d a;
0920 double b;
0921
0922 WorldToPolarCnv(s,polar,&a);
0923 d->az=a.x;
0924 d->el=a.y;
0925 d->r=a.z;
0926 b=tan(a.y);
0927 d->h=b*sin(a.x); // right plus
0928 d->v=b*cos(a.x); // up plus
0929}
0930
0931/**************************************************************************
0932 ClipAngleMatrix
0933
0934 Represent a side of a clip rect by a polar coordinate
0935 and an angular segment in its equatorial plane.
0936
0937 Input parameters are in the polar coordinate, where z axis
0938 is directed to the boresight.
0939
0940 azi azimuth (rad)
0941 inc inclination (rad - measured from z axis)
0942
0943 m conversion matrix (AngleMatrix)
0944 ang azimuth 0 to ang is a clip arc
0945**************************************************************************/
0946void ClipAngleMatrix(double azi,double inc,AngleMatrix *m,double *ang)
0947{
0948 double rot,az,de;
0949 double pi2=1.5707963268;
0950 Point2d s;
0951 Point3d d;
0952
0953 de=pi2-inc; // declination
0954
0955 m->c[0][0]=cos(de)*cos(azi); // body x' axis
0956 m->c[1][0]=cos(de)*sin(azi);
0957 m->c[2][0]=sin(de);
0958
0959 // rotate azimuth by an angle rot so that az falls 0 and 1.57079
0960 rot=0.;
0961 az=azi;
0962 while(az>pi2){
0963 az-=pi2;
0964 rot-=pi2;
0965 }
0966 while(az<0.){
0967 az+=pi2;
0968 rot+=pi2;
0969 }
0970
0971 // obtain *ang as angle between 2 vectors with the same declination
0972 d.x=cos(de)*cos(az);
0973 d.y=cos(de)*sin(az);
0974 d.z=sin(de);
0975 *ang=acos(d.x*d.x-d.y*d.y+d.z*d.z); // d and another vector symmet with x-xis
0976
0977 // now resume calculating z' axis
0978 s.x=-rot; // azimuth -rot
0979 s.y=-atan(tan(inc)*cos(az)); // elevation
0980 ToCartesian3d(&s,&d);
0981
0982 m->c[0][2]=d.x;
0983 m->c[1][2]=d.y;
0984 m->c[2][2]=d.z;
0985
0986 // y' axis as a vector product of z' and x'
0987 m->c[0][1]=m->c[1][2]*m->c[2][0] -m->c[2][2]*m->c[1][0];
0988 m->c[1][1]=m->c[2][2]*m->c[0][0] -m->c[0][2]*m->c[2][0];
0989 m->c[2][1]=m->c[0][2]*m->c[1][0] -m->c[1][2]*m->c[0][0];
0990
0991
0992}
0993
0994
0995/**************************************************************************
0996 SurfthePolarTAM
0997
0998 Make a coordinate system which converts two dimensional surface points
0999 to two dimensional polar coordinates on the viewer's sphere
1000
1001 s surface coordinate
1002 p polar coordinate ( usually 'thePolar' )
1003
1004 sp TransAngleMatrix where the surface is the body,
1005 and the polar is the frame
1006**************************************************************************/
1007void SurfthePolarTAM(TransAngleEuler *s,TransAngleMatrix *p,TransAngleMatrix *sp)
1008{
1009 TransAngleMatrix m;
1010 TransAngleEuler e,d;
1011
1012 InverseTAM(p,&m); // inverse of thePolar
1013 TAMtoTAE(&m,&e);
1014 MoveTwiceTAE(&e,s,&d);
1015 TAEtoTAM(&d,sp);
1016}
1017
1018/**************************************************************************
1019 SurfToPolar
1020
1021 Converts a two dimensional surface point
1022 to a two dimensional polar coordinate on the viewer's sphere
1023
1024 s surface coordinate (x,y)
1025 sp TransAngleMatrix where the surface is the body,
1026 and the polar is the frame
1027
1028 p polar coordinate (alpha,delta)
1029**************************************************************************/
1030void SurfToPolar(Point2d *s,TransAngleMatrix *sp,Point2d *p)
1031{
1032 Point3d a;
1033 double d;
1034
1035 a.x=sp->x+sp->c[0][0]*s->x+sp->c[0][1]*s->y;
1036 a.y=sp->y+sp->c[1][0]*s->x+sp->c[1][1]*s->y;
1037 a.z=sp->z+sp->c[2][0]*s->x+sp->c[2][1]*s->y;
1038
1039 d=sqrt(a.x*a.x+a.y*a.y);
1040 p->y=atan2(a.z,d);
1041 if(d==0.){
1042 p->x=0.;
1043 return;
1044 }
1045 p->x=atan2(a.y,a.x);
1046}
1047
1048/**************************************************************************
1049 PolarToSurf
1050
1051 Converts a two dimensional polar coordinate on the viewer's sphere
1052 to a two dimensional surface point
1053
1054 p polar coordinate (alpha,delta)
1055 sp TransAngleMatrix where the surface is the body,
1056 and the polar is the frame
1057
1058 s surface coordinate (x,y)
1059**************************************************************************/
1060void PolarToSurf(Point2d *p,TransAngleMatrix *sp,Point2d *s)
1061{
1062 Point3d a,b,c;
1063 double d;
1064
1065 a.x=cos(p->y)*cos(p->x);
1066 a.y=cos(p->y)*sin(p->x);
1067 a.z=sin(p->y);
1068 FrameToBodyTAM(&a,sp,&b);
1069 c.x=0.;
1070 c.y=0.;
1071 c.z=0.;
1072 FrameToBodyTAM(&c,sp,&a);
1073 s->x=(b.z*a.x-a.z*b.x)/(b.z-a.z);
1074 s->y=(b.z*a.y-a.z*b.y)/(b.z-a.z);
1075}
1076
1077/**************************************************************************
1078 CrossClipSide
1079
1080 Tests whether an arc given by 2 points in plolar coordinate
1081 intersects a side of the clipping rect, returnes if the arc
1082 lies outside of the clipping arc, if the arc could be inside
1083 of the clipping arc, or if the arc intersects the clipping arc.
1084 In the lattermost case, polar coordinate of the intersecting
1085 point will be returned as the last parameter.
1086
1087 frame coordinate system where a side of the clipping rect
1088 lies on the x-y plane, and x-axis corresponds to the
1089 start point of the side
1090 arc angle between the starting point and the end point
1091 of the arc corresponding to the clipping side
1092 p1 starting point the arc to be tested
1093 p2 end point the arc to be tested
1094 p crossing point
1095
1096 returned int 0 - could be inside of the clipping rect
1097 1 - intersect (p1 out, p2 in)
1098 2 - intersect (p1 in, p2 out)
1099 3 - outside
1100
1101**************************************************************************/
1102int CrossClipSide(AngleMatrix *frame,double arc,Point2d *p1,
1103 Point2d *p2,Point2d *p,double *ang)
1104{
1105 Point2d p3,p4,p8;
1106 Point3d p5,p6,p7;
1107
1108 FrameToBodyAMP(p1,frame,&p3);
1109 FrameToBodyAMP(p2,frame,&p4);
1110 if(p3.y>0.&&p4.y>0.) return 3; // outside
1111 if(p3.y<0.&&p4.y<0.) return 0;
1112 ToCartesian3d(&p3,&p5);
1113 ToCartesian3d(&p4,&p6);
1114
1115 VectorProduct(&p5,&p6,&p7);
1116
1117 ToPolar3d(&p7,&p8);
1118 if(p3.y>0.) p8.x-=1.5707963268; // p1 out, p2 in
1119 else p8.x+=1.5707963268;
1120
1121 if(p8.x>6.2831853072) p8.x-=6.2831853072;
1122 if(p8.x<0.) p8.x+=6.2831853072;
1123
1124 if(p8.x>arc) return 0; // could intersect with another side
1125
1126 p8.y=0.;
1127 BodyToFrameAMP(&p8,frame,p);
1128 *ang=p8.x; // angle from the corner
1129 if(p3.y>0.) return 1;
1130 else return 2;
1131}
1132
1133/****************************************************************************
1134 compactFarNear
1135
1136****************************************************************************/
1137/*
1138 void compactFarNear(StartPTE *farNear,int *n)
1139 {
1140 StartPTE a;
1141 int i,j;
1142 bool recursFarNear(StartPTE *farNear,int n,int p1,int p2,int depth);
1143
1144 for(i=0;i<*n;i++){
1145 a=farNear[i];
1146 farNear[i].start=-1;
1147 farNear[i].pte=-1;
1148 if(recursFarNear(farNear,*n,a.start,a.pte,0)){
1149 for(j=i+1;j<*n;j++) farNear[j-1]=farNear[j];
1150 *n--;
1151 }
1152 else farNear[i]=a;
1153 }
1154 }
1155*/
1156
1157
1158
1159/****************************************************************************
1160 inferFarNear
1161
1162 Given a set of far-near relationships, infers the far-near relationship
1163 of a given pair of items
1164
1165 farNear array of far-near relationhips - start:far, pte:near
1166 n number of data in array farNear
1167 p1 the first item
1168 p2 the second item
1169
1170 *near either p1 or p2, which is nearer
1171 returned if false, the relationship cannot be inferred
1172*****************************************************************************/
1173bool inferFarNear(StartPTE *farNear,long n,int p1,int p2,int *near)
1174{
1175
1176 // test if p1 is farther
1177 *near=p2;
1178 if(recursFarNear(farNear,n,p1,p2,0)) return true;
1179 // test if p2 is farther
1180 *near=p1;
1181 if(recursFarNear(farNear,n,p2,p1,0)) return true;
1182 return false;
1183}
1184
1185/****************************************************************************
1186 recursFarNear
1187
1188 finds recursively a path to the destination
1189
1190 farNear array of far-near relationhips - start:far, pte:near
1191 n number of data in array farNear
1192 p1 the first item
1193 p2 the second item
1194 depth depth of recursion - detect a loop condition and abort
1195
1196 returned if false, the relationship cannot be inferred
1197*****************************************************************************/
1198bool recursFarNear(StartPTE *farNear,long n,int p1,int p2,long depth)
1199{
1200 long i;
1201
1202 depth++;
1203 if(depth>n) return false; // escape loop
1204
1205 for(i=0;i<n;i++){
1206 if(farNear[i].start==p1){
1207 if(farNear[i].pte==p2) return true;
1208 if(recursFarNear(farNear,n,farNear[i].pte,p2,depth)) return true;
1209 }
1210 }
1211 return false;
1212}
1213
1214/****************************************************************************
1215 Overlap
1216
1217 tests if the given two surfaces are seen overlapped
1218
1219 p point array
1220 s1 indeces of points for a surface 1
1221 s2 indices of points for a surface 2
1222
1223 p0 typical overlapped point
1224
1225 returned if false, the surfaces are disjoint
1226*****************************************************************************/
1227bool Overlap(Point2d *p,StartPTE *s1,StartPTE *s2,Point2d *p0)
1228{
1229 Point2d s1max,s1min;
1230 Point2d s2max,s2min;
1231 int i,j;
1232 double x1,y1,x2,y2,x3,y3,x4,y4;
1233 double a,b,d;
1234 double l1,l2;
1235
1236 // for(i=s1->start;i<s1->pte;i++){
1237 // sprintf(msg,"Overlap s1 i= %3i p[i].x=%8.3f p[i].y=%8.3f\n",i,p[i].x,p[i].y);
1238 // outstr(msg);
1239 //}
1240 //flush();
1241 //for(i=s2->start;i<s2->pte;i++){
1242 // sprintf(msg,"Overlap s2 i= %3i p[i].x=%8.3f p[i].y=%8.3f\n",i,p[i].x,p[i].y);
1243 // outstr(msg);
1244 //}
1245 //flush();
1246
1247 s1max.x=-10000.;
1248 s1max.y=-10000.;
1249 s1min.x=10000.;
1250 s1min.y=10000.;
1251
1252 for(i=s1->start;i<s1->pte;i++){
1253 if(p[i].x>s1max.x) s1max.x=p[i].x;
1254 if(p[i].y>s1max.y) s1max.y=p[i].y;
1255 if(p[i].x<s1min.x) s1min.x=p[i].x;
1256 if(p[i].y<s1min.y) s1min.y=p[i].y;
1257 }
1258
1259 s2max.x=-10000.;
1260 s2max.y=-10000.;
1261 s2min.x=10000.;
1262 s2min.y=10000.;
1263
1264 for(i=s2->start;i<s2->pte;i++){
1265 if(p[i].x>s2max.x) s2max.x=p[i].x;
1266 if(p[i].y>s2max.y) s2max.y=p[i].y;
1267 if(p[i].x<s2min.x) s2min.x=p[i].x;
1268 if(p[i].y<s2min.y) s2min.y=p[i].y;
1269 }
1270
1271 if(s1max.x<s2min.x) return false;
1272 if(s2max.x<s1min.x) return false;
1273 if(s1max.y<s2min.y) return false;
1274 if(s2max.y<s1min.y) return false;
1275
1276 for(i=s1->start;i<s1->pte;i++){
1277 x1=p[i].x;
1278 y1=p[i].y;
1279 if(i==s1->pte-1){
1280 if(s1->pte-s1->start==2) break;
1281 x2=p[s1->start].x;
1282 y2=p[s1->start].y;
1283 }
1284 else{
1285 x2=p[i+1].x;
1286 y2=p[i+1].y;
1287 }
1288 l1=(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1); // approximate l1*l2 by this
1289 for(j=s2->start;j<s2->pte;j++){
1290 x3=p[j].x;
1291 y3=p[j].y;
1292 if(j==s2->pte-1){
1293 if(s2->pte-s2->start==2) break;
1294 x4=p[s2->start].x;
1295 y4=p[s2->start].y;
1296 }
1297 else{
1298 x4=p[j+1].x;
1299 y4=p[j+1].y;
1300 }
1301 // detect simply disjoint cases
1302 while(1){
1303 if(x1<=x2){
1304 if(x3<=x4){
1305 if(x2<x3) break;
1306 if(x4<x1) break;
1307 }
1308 else{
1309 if(x2<x4) break;
1310 if(x3<x1) break;
1311 }
1312 }
1313 else{
1314 if(x3<=x4){
1315 if(x1<x3) break;
1316 if(x4<x2) break;
1317 }
1318 else{
1319 if(x1<x4) break;
1320 if(x3<x2) break;
1321 }
1322 }
1323 if(y1<=y2){
1324 if(y3<=y4){
1325 if(y2<y3) break;
1326 if(y4<y1) break;
1327 }
1328 else{
1329 if(y2<y4) break;
1330 if(y3<y1) break;
1331 }
1332 }
1333 else{
1334 if(y3<=y4){
1335 if(y1<y3) break;
1336 if(y4<y2) break;
1337 }
1338 else{
1339 if(y1<y4) break;
1340 if(y3<y2) break;
1341 }
1342 }
1343 // sprintf(msg,"Overlap x1=%8.3f y1=%8.3f x2=%8.3f y2=%8.3f x3=%8.3f y3=%8.3f x4=%8.3f y4=%8.3f \n",x1,y1,x2,y2,x3,y3,x4,y4);
1344 // outstr(msg);
1345 // flush();
1346 d=(x2-x1)*(y3-y4)-(x3-x4)*(y2-y1);
1347 if(fabs(d/l1)<0.001) break;
1348 a=((x3-x1)*(y3-y4)-(x3-x4)*(y3-y1))/d;
1349 // if(a<-0.01||a>1.01) break;
1350 if(a<0.01||a>0.99) break;
1351 b=((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1))/d;
1352 // if(b<-0.01||b>1.01) break;
1353 if(b<0.01||b>0.99) break;
1354 if(a<0.01&&b<0.01) break;
1355 if(a<0.01&&b>0.99) break;
1356 if(a>0.99&&b<0.01) break;
1357 if(a>0.99&&b>0.99) break;
1358 p0->x=x1+(x2-x1)*a;
1359 p0->y=y1+(y2-y1)*a;
1360 return true;
1361 }
1362 }
1363 }
1364 // test which bounding rect encloses the other bounding rect
1365 if((s1max.x-s1min.x)<=(s2max.x-s2min.x)&&
1366 (s1max.y-s1min.y)<=(s2max.y-s2min.y)){
1367 i=s1->start;
1368 x1=0.5*(p[i].x+p[i+2].x);
1369 y1=0.5*(p[i].y+p[i+2].y);
1370 d=0.;
1371 for(j=s2->start;j<s2->pte;j++){
1372 x3=p[j].x;
1373 y3=p[j].y;
1374 if(j==s2->pte-1){
1375 if(s2->pte-s2->start==2) break;
1376 x4=p[s2->start].x;
1377 y4=p[s2->start].y;
1378 }
1379 else{
1380 x4=p[j+1].x;
1381 y4=p[j+1].y;
1382 }
1383 d+=ArgDif(x1,y1,x3,y3,x4,y4);
1384 }
1385 if(d>6.){
1386 p0->x=x1;
1387 p0->y=y1;
1388 return true;
1389 }
1390
1391 }
1392
1393 if((s1max.x-s1min.x)>=(s2max.x-s2min.x)&&
1394 (s1max.y-s1min.y)>=(s2max.y-s2min.y)){
1395 i=s2->start;
1396 x1=0.5*(p[i].x+p[i+2].x);
1397 y1=0.5*(p[i].y+p[i+2].y);
1398 d=0.;
1399 for(j=s1->start;j<s1->pte;j++){
1400 x3=p[j].x;
1401 y3=p[j].y;
1402 if(j==s1->pte-1){
1403 if(s1->pte-s1->start==2) break;
1404 x4=p[s1->start].x;
1405 y4=p[s1->start].y;
1406 }
1407 else{
1408 x4=p[j+1].x;
1409 y4=p[j+1].y;
1410 }
1411 d+=ArgDif(x1,y1,x3,y3,x4,y4);
1412 }
1413 if(d>6.){
1414 p0->x=x1;
1415 p0->y=y1;
1416 return true;
1417 }
1418 }
1419 return false;
1420}
1421
1422/****************************************************************************
1423 FarNearDist
1424
1425 distance separating surface1 and surface2 for a given direction of view
1426
1427 s1 coordinate conversion (surface:body view:frame)
1428 s2 coordinate conversion (surface:body view:frame)
1429 p polar coordinate of the line of sight
1430
1431 d distance to s1 minus distance to s2
1432*****************************************************************************/
1433void FarNearDist(TransAngleMatrix *s1,TransAngleMatrix *s2,Point2d *p,double *d)
1434{
1435 Point2d srf1,srf2;
1436 Point3d a;
1437
1438 PolarToSurf(p,s1,&srf1);
1439 PolarToSurf(p,s2,&srf2);
1440
1441 a.x=s1->x+s1->c[0][0]*srf1.x+s1->c[0][1]*srf1.y;
1442 a.y=s1->y+s1->c[1][0]*srf1.x+s1->c[1][1]*srf1.y;
1443 a.z=s1->z+s1->c[2][0]*srf1.x+s1->c[2][1]*srf1.y;
1444
1445 *d=sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
1446
1447 a.x=s2->x+s2->c[0][0]*srf2.x+s2->c[0][1]*srf2.y;
1448 a.y=s2->y+s2->c[1][0]*srf2.x+s2->c[1][1]*srf2.y;
1449 a.z=s2->z+s2->c[2][0]*srf2.x+s2->c[2][1]*srf2.y;
1450
1451 *d-=sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
1452
1453}
1454
1455/****************************************************************************
1456 PolarToGnomo
1457
1458 Polar to Gnomonic Conversion
1459
1460 p polar coordinate (x:aximuth, y:elevation)
1461 g gnomonic coordinate (x:toward polar y, y:toward polar x)
1462*****************************************************************************/
1463void PolarToGnomo(Point2d *p,Point2d *g)
1464{
1465 double r;
1466
1467 r=tan(1.5707963267-p->y);
1468 g->x=r*sin(p->x);
1469 g->y=r*cos(p->x);
1470}
1471
1472
1473/****************************************************************************
1474 GnomoToPolar
1475
1476 Gnomonic to Polar Conversion
1477
1478 g gnomonic coordinate (x:toward polar y, y:toward polar x)
1479 p polar coordinate (x:aximuth, y:elevation)
1480*****************************************************************************/
1481void GnomoToPolar(Point2d *g,Point2d *p)
1482{
1483 double r;
1484
1485 r=sqrt(g->x*g->x+g->y*g->y);
1486 p->y=1.5707963267-atan(r);
1487 p->x=atan2(g->x,g->y);
1488}
1489
1490
1491/****************************************************************************
1492 ArgDif
1493
1494 argument difference between two vectors
1495
1496 (x0,y0) origin
1497 (x1,y1) vector1 tip
1498 (x2,y2) vector2 tip
1499
1500 returned argument of vector2 - argument of vector1
1501*****************************************************************************/
1502double ArgDif(double x0,double y0,double x1,double y1,double x2,double y2)
1503{
1504 double ang;
1505
1506 ang=atan2(y2-y0,x2-x0)-atan2(y1-y0,x1-x0);
1507 if(fabs(ang)<3.141592) return ang;
1508 if(ang<0.) ang+=6.283185;
1509 else ang-=6.283185;
1510 return ang;
1511}
1512
1513
1514/****************************************************************************
1515 SameSideOfSurface
1516
1517 Return true if the Two Points are in the Same Side of a Surface
1518
1519 s Surface
1520 p1 point 1
1521 p2 point 2
1522
1523 returned true if point 1 and 2 are in the same side of s
1524*****************************************************************************/
1525bool SameSideOfSurface(Surface *s,Point3d *p1,Point3d *p2)
1526{
1527 Point3d z,d1,d2;
1528
1529 z.x=sin(s->theta0)*sin(s->phi0); // axis normal to the plane
1530 z.y=-sin(s->theta0)*cos(s->phi0);
1531 z.z=cos(s->theta0);
1532 d1.x=p1->x-s->x0;
1533 d1.y=p1->y-s->y0;
1534 d1.z=p1->z-s->z0;
1535 d2.x=p2->x-s->x0;
1536 d2.y=p2->y-s->y0;
1537 d2.z=p2->z-s->z0;
1538 if((d1.x*z.x+d1.y*z.y+d1.z*z.z)*(d2.x*z.x+d2.y*z.y+d2.z*z.z)>0.)
1539 return true;
1540 else return false;
1541}
1542
1543
1544
1545/******************************************************************************
1546 CMPpart
1547
1548 part of the concept
1549******************************************************************************/
1550void CMPpart::setCood(double x,double y,double z,double theta,double phi,double psi)
1551{
1552 x0=x; /* local coordinate referred to */
1553 y0=y; /* refPart */
1554 z0=z;
1555 theta0=theta;
1556 phi0=phi;
1557 psi0=psi;
1558}
1559
1560void CMPpart::setColor(int i)
1561{
1562 color=i;
1563}
1564
1565void CMPpart::hide()
1566{
1567 showDrawing=false;
1568}
1569
1570
1571/******************************************************************************
1572 CMPpoint
1573
1574 part of the concept - point
1575******************************************************************************/
1576
1577
1578/******************************************************************************
1579 CMPlines
1580
1581 part of the concept - lines
1582******************************************************************************/
1583
1584
1585void CMPlines::SetPoint(double x,double y,double z)
1586{
1587 Point3d p;
1588 p.x=x;
1589 p.y=y;
1590 p.z=z;
1591 pt.push_back(p);
1592 npoint++;
1593}
1594
1595
1596/******************************************************************************
1597 GetSurface
1598
1599 Return a line segment
1600
1601******************************************************************************/
1602
1603void CMPlines::GetSurface(int m,TransAngleMatrix *po,Surface *surf,Point2d *surfPt)
1604{
1605 double x1,y1,z1,x2,y2,z2;
1606 TransAngleEuler tae1;
1607 TransAngleMatrix tam1;
1608 TransAngleEuler tae2;
1609 TransAngleEuler tae;
1610 Point3d s;
1611 Point3d e;
1612 Point3d f;
1613 Point3d g;
1614 double length,d;
1615 Point2d a;
1616
1617 s=pt[m];
1618 e.x=pt[m+1].x-s.x;
1619 e.y=pt[m+1].y-s.y;
1620 e.z=pt[m+1].z-s.z;
1621
1622 if(e.x==0.&&e.y==0.&&e.z==0.){
1623 surf->start=-1;
1624 surf->pte=-1;
1625 return;
1626 }
1627
1628 tae1.x=xg; // Translate Angle Euler for the part
1629 tae1.y=yg;
1630 tae1.z=zg;
1631 tae1.theta=thetag;
1632 tae1.phi=phig;
1633 tae1.psi=psig;
1634 TAEtoTAM(&tae1,&tam1);
1635 FrameToBodyTAM((Point3d *)po,&tam1,&g); // viewer by part body coordinate
1636 g.x-=s.x; // direction to the viewer
1637 g.y-=s.y;
1638 g.z-=s.z;
1639 NormPoint3d(&g,&g,&d);
1640 tam1.x=s.x; // Translate Angle Matrix for the line segment
1641 tam1.y=s.y;
1642 tam1.z=s.z;
1643 NormPoint3d(&e,&f,&length); // normalized y-axis
1644 if(fabs(f.x*g.x+f.y*g.y+f.z*g.z)>0.99996){ // if the line is directed to viewer
1645 ToPolar3d(&g,&a); // viewer's direction
1646 a.y-=0.0174532; // Offset 1 deg from the direction to the viewer
1647 f.x=cos(a.y)*cos(a.x); // updated line direction to avoid anomaly
1648 f.y=cos(a.y)*sin(a.x);
1649 f.z=sin(a.y);
1650 }
1651 tam1.c[0][1]=f.x;
1652 tam1.c[1][1]=f.y;
1653 tam1.c[2][1]=f.z;
1654 VectorProduct(&f,&g,&s); // direction to x-axis
1655 NormPoint3d(&s,&s,&d);
1656 tam1.c[0][0]=s.x;
1657 tam1.c[1][0]=s.y;
1658 tam1.c[2][0]=s.z;
1659 VectorProduct(&s,&f,&g); // direction to z-axis
1660 tam1.c[0][2]=g.x;
1661 tam1.c[1][2]=g.y;
1662 tam1.c[2][2]=g.z;
1663 TAMtoTAE(&tam1,&tae2);
1664 MoveTwiceTAE(&tae1,&tae2,&tae);
1665 surf->x0=tae.x;
1666 surf->y0=tae.y;
1667 surf->z0=tae.z;
1668 surf->theta0=tae.theta;
1669 surf->phi0=tae.phi;
1670 surf->psi0=tae.psi;
1671 surf->start=0;
1672 surf->pte=2;
1673 surfPt[0].x=0.;
1674 surfPt[0].y=0.;
1675 surfPt[1].x=0.;
1676 surfPt[1].y=length;
1677}
1678
1679/******************************************************************************
1680 GetDrawData
1681
1682 Return data for drawing
1683
1684 m line segment number (0-)
1685 po thePolar
1686
1687 seg array of line segments (clockwise angle from top, angle from line of sight)
1688 blk indices of segs which are to be drawn with lines
1689 cont indices of closed line segments. Inside of the closure
1690 will be painted by a color.
1691
1692******************************************************************************/
1693
1694void CMPlines::GetDrawData(int m,TransAngleMatrix *po,
1695 Line2d *seg,StartPTE *blk,StartPTE *cont)
1696{
1697 TransAngleEuler tae;
1698 TransAngleMatrix tam;
1699 Point3d s;
1700 Point3d e;
1701 Point3d f;
1702 double r;
1703 Point2d a;
1704
1705 s=pt[m];
1706 e=pt[m+1];
1707
1708 tae.x=xg; // Translate Angle Euler for the part
1709 tae.y=yg;
1710 tae.z=zg;
1711 tae.theta=thetag;
1712 tae.phi=phig;
1713 tae.psi=psig;
1714 TAEtoTAM(&tae,&tam);
1715 BodyToFrameTAM(&s,&tam,&f);
1716 FrameToBodyTAM(&f,po,&s);
1717 NormPoint3d(&s,&f,&r);
1718 ToPolar3d(&f,&a);
1719 a.y=1.5707963267-a.y;
1720 seg[0].sx=a.x;
1721 seg[0].sy=a.y;
1722 BodyToFrameTAM(&e,&tam,&f);
1723 FrameToBodyTAM(&f,po,&e);
1724 NormPoint3d(&e,&f,&r);
1725 ToPolar3d(&f,&a);
1726 a.y=1.5707963267-a.y;
1727 seg[0].ex=a.x;
1728 seg[0].ey=a.y;
1729 blk->start=0;
1730 blk->pte=1;
1731 cont->start=0;
1732 cont->pte=1;
1733
1734}
1735
1736
1737/******************************************************************************
1738 FieldOfView
1739
1740 Get horizontal and vertical field of view to cover the object
1741******************************************************************************/
1742
1743void CMPlines::FieldOfView(TransAngleMatrix *vu,Point2d *az,Point2d *el)
1744{
1745 TransAngleEuler tae;
1746 TransAngleMatrix tam;
1747 int i;
1748 Point3d p;
1749 Point3d s;
1750 Point3d f;
1751 Point2d a;
1752 double d;
1753
1754 tae.x=xg; // Translate Angle Euler for the part
1755 tae.y=yg;
1756 tae.z=zg;
1757 tae.theta=thetag;
1758 tae.phi=phig;
1759 tae.psi=psig;
1760 TAEtoTAM(&tae,&tam);
1761 az->x=10.;
1762 az->y=-10.;
1763 el->x=10.;
1764 el->y=-10.;
1765 for(i=0;i<npoint;i++){
1766 s=pt[i];
1767 BodyToFrameTAM(&s,&tam,&f);
1768 FrameToBodyTAM(&f,vu,&s);
1769 NormPoint3d(&s,&s,&d);
1770 ToPolar3d(&s,&a);
1771 if(a.x<az->x) az->x=a.x;
1772 if(a.x>az->y) az->y=a.x;
1773 if(a.y<el->x) el->x=a.y;
1774 if(a.y>el->y) el->y=a.y;
1775 }
1776}
1777
1778/******************************************************************************
1779 CMPplane
1780
1781 part of the concept - plane
1782******************************************************************************/
1783
1784void CMPplane::SetPoint(double x,double y)
1785{
1786 Point2d p;
1787 p.x=x;
1788 p.y=y;
1789 pt.push_back(p);
1790 npoint++;
1791}
1792
1793/******************************************************************************
1794 GetSurface
1795
1796 Return the surface data
1797
1798 The surface is assumed to be given by points that are arranged to
1799 move clock-wise looking toward the z vector.
1800
1801 Both sides of the surface are assumed to be visible. Therefore,
1802 the order of the point array is reversed, when the viewer's eye
1803 is in the -z side.
1804
1805******************************************************************************/
1806
1807void CMPplane::GetSurface(TransAngleMatrix *po,Surface *surf,Point2d *surfPt)
1808{
1809 TransAngleMatrix tam;
1810 Point3d d;
1811 int i,j;
1812 int np;
1813 double x1,y1;
1814
1815 surf->x0=xg;
1816 surf->y0=yg;
1817 surf->z0=zg;
1818 surf->theta0=thetag;
1819 surf->phi0=phig;
1820 surf->psi0=psig;
1821 surf->start=0;
1822 TAEtoTAM((TransAngleEuler *)surf,&tam);
1823 FrameToBodyTAM((Point3d *)po,&tam,&d); // if d>0 +z axis is facing the viewer
1824 np=0; // number of effective points
1825 if(d.z>0.){
1826 x1=pt[npoint-1].x;
1827 y1=pt[npoint-1].y;
1828 }
1829 else{
1830 x1=pt[0].x;
1831 y1=pt[0].y;
1832 }
1833 for(i=0;i<npoint;i++){
1834 if(d.z>0.) j=i;
1835 else j=npoint-i-1;
1836 surfPt[np].x=pt[j].x;
1837 surfPt[np].y=pt[j].y;
1838 if(x1!=surfPt[np].x||y1!=surfPt[np].y){
1839 x1=surfPt[np].x;
1840 y1=surfPt[np].y;
1841 if(np<47) np++; // taking into some allowance to max 50
1842 }
1843 }
1844 surf->pte=np;
1845}
1846
1847
1848/******************************************************************************
1849 GetDrawData
1850
1851 Return data for drawing
1852
1853 po thePolar
1854
1855 seg array of line segments (clockwise angle from top, angle from line of sight)
1856 blk indices of segs which are to be drawn with lines
1857 cont indices of closed line segments. Inside of the closure
1858 will be painted by a color.
1859
1860******************************************************************************/
1861
1862void CMPplane::GetDrawData(TransAngleMatrix *po,Line2d *seg,StartPTE *blk,StartPTE *cont)
1863{
1864 StartPTE blk1;
1865 StartPTE cont1;
1866 int i,j;
1867 double sx1,sy1,ex1,ey1;
1868 CMPplane *next;
1869
1870 cont->start=0;
1871 GetDrawData1(po,seg,blk,cont);
1872 if(otherHalf){
1873 bool face1,face2;
1874 face1=FacingViewer((Point3d *)po);
1875 next=otherHalf;
1876 while(next!=this){
1877 cont1.start=blk->pte;
1878 next->GetDrawData1(po,seg,&blk1,&cont1);
1879 face2=next->FacingViewer((Point3d *)po);
1880 for(i=blk->start;i<cont->pte;i++){
1881 if(face1==face2){
1882 sx1=seg[i].sx;
1883 sy1=seg[i].sy;
1884 ex1=seg[i].ex;
1885 ey1=seg[i].ey;
1886 }
1887 else{
1888 sx1=seg[i].ex;
1889 sy1=seg[i].ey;
1890 ex1=seg[i].sx;
1891 ey1=seg[i].sy;
1892 }
1893 for(j=cont1.start;j<cont1.pte;j++){
1894 double d0,d1,d2,d3;
1895 d0=seg[j].sx-ex1;
1896 d1=seg[j].sy-ey1;
1897 d2=seg[j].ex-sx1;
1898 d3=seg[j].ey-sy1;
1899 if(d0*d0+d1*d1+d2*d2+d3*d3<0.000001) break;
1900 // if(seg[j].sx==ex1&&seg[j].sy==ey1&&seg[j].ex==sx1&&seg[j].ey==sy1)
1901 // break;
1902 }
1903 if(j!=cont1.pte){
1904 for(j=i;j>0;j--) seg[j]=seg[j-1];
1905 seg[0].sx=sx1;
1906 seg[0].sy=sy1;
1907 seg[0].ex=ex1;
1908 seg[0].ey=ey1;
1909 blk->start++;
1910 }
1911 }
1912 next=next->otherHalf;
1913 }
1914 }
1915
1916
1917}
1918
1919void CMPplane::GetDrawData1(TransAngleMatrix *po,Line2d *seg,StartPTE *blk,StartPTE *cont)
1920{
1921 TransAngleEuler tae;
1922 TransAngleMatrix tam;
1923 Point3d a,b,c,d;
1924 Point2d e,f,g;
1925 double r;
1926 int i,j,k;
1927
1928 k=cont->start;
1929 tae.x=xg; // Translate Angle Euler for the part
1930 tae.y=yg;
1931 tae.z=zg;
1932 tae.theta=thetag;
1933 tae.phi=phig;
1934 tae.psi=psig;
1935 TAEtoTAM(&tae,&tam);
1936 a.x=pt[0].x;
1937 a.y=pt[0].y;
1938 a.z=0.;
1939 BodyToFrameTAM(&a,&tam,&b);
1940 FrameToBodyTAM(&b,po,&c);
1941 NormPoint3d(&c,&d,&r);
1942 ToPolar3d(&d,&e);
1943 e.y=1.5707963267-e.y;
1944 f=e;
1945 for(i=1;i<=npoint;i++){
1946 g=f;
1947 if(i<npoint){
1948 a.x=pt[i].x;
1949 a.y=pt[i].y;
1950 a.z=0.;
1951 BodyToFrameTAM(&a,&tam,&b);
1952 FrameToBodyTAM(&b,po,&c);
1953 NormPoint3d(&c,&d,&r);
1954 ToPolar3d(&d,&f);
1955 f.y=1.5707963267-f.y;
1956 }
1957 else f=e;
1958 seg[k].sx=g.x;
1959 seg[k].sy=g.y;
1960 seg[k].ex=f.x;
1961 seg[k].ey=f.y;
1962 k++;
1963
1964 }
1965 blk->start=cont->start;
1966 blk->pte=blk->start+npoint;
1967 cont->pte=cont->start+npoint;
1968}
1969
1970/******************************************************************************
1971 FieldOfView
1972
1973 Get horizontal and vertical field of view to cover the object
1974******************************************************************************/
1975
1976void CMPplane::FieldOfView(TransAngleMatrix *vu,Point2d *az,Point2d *el)
1977{
1978 TransAngleEuler tae;
1979 TransAngleMatrix tam;
1980 int i;
1981 Point2d p;
1982 Point3d s;
1983 Point3d f;
1984 Point2d a;
1985 double d;
1986
1987 tae.x=xg; // Translate Angle Euler for the part
1988 tae.y=yg;
1989 tae.z=zg;
1990 tae.theta=thetag;
1991 tae.phi=phig;
1992 tae.psi=psig;
1993 TAEtoTAM(&tae,&tam);
1994 az->x=10.;
1995 az->y=-10.;
1996 el->x=10.;
1997 el->y=-10.;
1998 for(i=0;i<npoint;i++){
1999 s.x=pt[i].x;
2000 s.y=pt[i].y;
2001 s.z=0.;
2002 BodyToFrameTAM(&s,&tam,&f);
2003 FrameToBodyTAM(&f,vu,&s);
2004 NormPoint3d(&s,&s,&d);
2005 ToPolar3d(&s,&a);
2006 if(a.x<az->x) az->x=a.x;
2007 if(a.x>az->y) az->y=a.x;
2008 if(a.y<el->x) el->x=a.y;
2009 if(a.y>el->y) el->y=a.y;
2010 }
2011}
2012
2013/*****************************************************************
2014 Return if the surface is facing the viewer
2015
2016 po thePolar casted to Point3d
2017*****************************************************************/
2018bool CMPplane::FacingViewer(Point3d *po)
2019{
2020 double nx,ny,nz;
2021
2022 ny=sin(thetag);
2023 nx=ny*cos(phig-1.5707963268);
2024 ny*=sin(phig-1.5707963268);
2025 nz=cos(thetag);
2026 return (po->x-xg)*nx+(po->y-yg)*ny+(po->z-zg)*nz>0.;
2027}
2028
2029
2030/******************************************************************************
2031 CMPbox
2032
2033******************************************************************************/
2034
2035void CMPbox::SetSize(double x,double y,double z)
2036{
2037 lx=x;
2038 ly=y;
2039 lz=z;
2040}
2041
2042
2043/******************************************************************************
2044 ClearVisSurface
2045
2046 Clear 'visibleSurface' instance variable so that new visibility
2047 from a new viewpoint can be obtained in GetSurface
2048******************************************************************************/
2049void CMPbox::ClearVisSurface(void)
2050{
2051 visibleSurface=0;
2052}
2053
2054/******************************************************************************
2055 GetSurface
2056
2057 Return a surface
2058
2059 Six surfaces are related to an index number in the following way
2060
2061 index face toward...
2062 0 +x
2063 1 +y
2064 2 +z
2065 3 -x
2066 4 -y
2067 5 -z
2068
2069
2070 The instance variable visibleSurface is used to map between
2071 the above index and the reference index passed as the first parameter.
2072
2073 The returned bool value is false if no more visible surface exists
2074 for the asked reference number.
2075******************************************************************************/
2076
2077bool CMPbox::GetSurface(int n,TransAngleMatrix *po,Surface *surf,Point2d *surfPt)
2078{
2079 double x1,y1,z1,x2,y2;
2080 TransAngleEuler tae1;
2081 TransAngleMatrix tam1;
2082 TransAngleEuler tae2;
2083 TransAngleEuler tae;
2084 Point3d p;
2085 Point3d s;
2086 Point3d e;
2087 Point3d f;
2088 double length;
2089 Point2d a;
2090 int i,j,m;
2091
2092 x1=0.5*lx;
2093 y1=0.5*ly;
2094 z1=0.5*lz;
2095 tae1.x=xg; // Translate Angle Euler for the part
2096 tae1.y=yg;
2097 tae1.z=zg;
2098 tae1.theta=thetag;
2099 tae1.phi=phig;
2100 tae1.psi=psig;
2101
2102 if(visibleSurface==0){
2103 TAEtoTAM(&tae1,&tam1);
2104 j=1;
2105 for(i=0;i<6;i++){
2106 s.x=((i>2)? -1.:1.)*((i%3==0)? x1:0.);
2107 s.y=((i>2)? -1.:1.)*((i%3==1)? y1:0.);
2108 s.z=((i>2)? -1.:1.)*((i%3==2)? z1:0.);
2109 BodyToFrameTAM(&s,&tam1,&e);
2110 if((xg-e.x)*(po->x-e.x)+(yg-e.y)*(po->y-e.y)
2111 +(zg-e.z)*(po->z-e.z)<0.) {
2112 // DebugStr("\psurface");
2113 visibleSurface|=j;
2114 }
2115 j=j<<1;
2116 }
2117 }
2118 j=1;
2119 i=0;
2120 for(m=0;m<6;m++){
2121 if(visibleSurface&j){
2122 if(i==n) break;
2123 i++;
2124 }
2125 j=j<<1;
2126 }
2127 if(m==6) return false;
2128 switch (m){
2129 case 0: // +x
2130 tae2.x=x1; // Translate Angle Euler for the side
2131 tae2.y=0.;
2132 tae2.z=0.;
2133 tae2.theta=1.5707963268; // rotate 90 deg around +y axis
2134 tae2.phi=1.5707963268;
2135 tae2.psi=-1.5707963268;
2136 x2=z1;
2137 y2=y1;
2138 break;
2139 case 1: // +y
2140 tae2.x=0.; // Translate Angle Euler for the side
2141 tae2.y=y1;
2142 tae2.z=0.;
2143 tae2.theta=1.5707963268; // rotate 90 deg around -x axis
2144 tae2.phi=3.141592653589;
2145 tae2.psi=3.141592653589;
2146 x2=x1;
2147 y2=z1;
2148 break;
2149 case 2: // +z
2150 tae2.x=0.; // Translate Angle Euler for the side
2151 tae2.y=0.;
2152 tae2.z=z1;
2153 tae2.theta=0.; // no rotation
2154 tae2.phi=0.;
2155 tae2.psi=0.;
2156 x2=x1;
2157 y2=y1;
2158 break;
2159 case 3: // -x
2160 tae2.x=-x1; // Translate Angle Euler for the side
2161 tae2.y=0.;
2162 tae2.z=0.;
2163 tae2.theta=1.5707963268; // rotate 90 deg around -y axis
2164 tae2.phi=-1.5707963268;
2165 tae2.psi=1.5707963268;
2166 x2=z1;
2167 y2=y1;
2168 break;
2169 case 4: // -y
2170 tae2.x=0.; // Translate Angle Euler for the side
2171 tae2.y=-y1;
2172 tae2.z=0.;
2173 tae2.theta=1.5707963268; // rotate 90 deg around -x axis
2174 tae2.phi=0.;
2175 tae2.psi=0.;
2176 x2=x1;
2177 y2=z1;
2178 break;
2179 case 5: // -z
2180 tae2.x=0.; // Translate Angle Euler for the side
2181 tae2.y=0.;
2182 tae2.z=-z1;
2183 tae2.theta=3.1415926536; // no rotation
2184 tae2.phi=0.;
2185 tae2.psi=0.;
2186 x2=x1;
2187 y2=y1;
2188 break;
2189 }
2190 MoveTwiceTAE(&tae1,&tae2,&tae);
2191 surf->x0=tae.x;
2192 surf->y0=tae.y;
2193 surf->z0=tae.z;
2194 surf->theta0=tae.theta;
2195 surf->phi0=tae.phi;
2196 surf->psi0=tae.psi;
2197 surf->start=0;
2198 surf->pte=4;
2199 surfPt[0].x=x2;
2200 surfPt[0].y=y2;
2201 surfPt[1].x=-x2;
2202 surfPt[1].y=y2;
2203 surfPt[2].x=-x2;
2204 surfPt[2].y=-y2;
2205 surfPt[3].x=x2;
2206 surfPt[3].y=-y2;
2207 return true;
2208}
2209
2210/******************************************************************************
2211 GetCenter
2212
2213 Return the center point in frame coordinate
2214
2215******************************************************************************/
2216
2217void CMPbox::GetCenter(Point3d *p)
2218{
2219 p->x=xg;
2220 p->y=yg;
2221 p->z=zg;
2222}
2223
2224/******************************************************************************
2225 GetDrawData
2226
2227 Return data for drawing
2228
2229 po thePolar
2230
2231 seg array of line segments (clockwise angle from top, angle from line of sight)
2232 blk indices of segs which are to be drawn with lines
2233 cont indices of closed line segments. Inside of the closure
2234 will be painted by a color.
2235
2236******************************************************************************/
2237
2238void CMPbox::GetDrawData(TransAngleMatrix *po,Line2d *seg,StartPTE *blk,StartPTE *cont)
2239{
2240 double x,y,z,x1,y1,z1,x2,y2,z2,dx,dy,dz;
2241 int isen[3][10];
2242 int i,j,k,l,m,n,nsen,n1,brnch,lin,line;
2243 TransAngleEuler tae1;
2244 TransAngleMatrix tam1;
2245 TransAngleEuler tae2;
2246 TransAngleEuler tae;
2247 Point3d apex;
2248 Point3d ten[8]; /* coordinate of 8 apexes */
2249 Point3d s,d;
2250 double r;
2251 Point2d direc;
2252
2253
2254 /* set conversion parameters */
2255 dx=lx;
2256 dy=ly;
2257 dz=lz;
2258
2259 tae1.x=xg; // Translate Angle Euler for the part
2260 tae1.y=yg;
2261 tae1.z=zg;
2262 tae1.theta=thetag;
2263 tae1.phi=phig;
2264 tae1.psi=psig;
2265 TAEtoTAM(&tae1,&tam1);
2266
2267 /* fill coordinate of the 8 apexes */
2268 /* l=0 ( 0, 0, 0) l=1 (dx, 0, 0) l=2 ( 0,dy, 0) l=3 (dx,dy, 0) */
2269 /* l=4 ( 0, 0,dz) l=5 (dx, 0,dz) l=6 ( 0,dy,dz) l=7 (dx,dy,dz) */
2270 for(i=0;i<2;i++){
2271 apex.z=dz*i-0.5*dz; /* float */
2272 for(j=0;j<2;j++){
2273 apex.y=dy*j-0.5*dy; /* float */
2274 for(k=0;k<2;k++){
2275 apex.x=dx*k-0.5*dx; /* float */
2276 l=4*i+2*j+k; /* l=0-7 */
2277 BodyToFrameTAM(&apex,&tam1,&ten[l]);
2278 }
2279 }
2280 }
2281
2282 /* fill line information */
2283 for(i=0;i<10;i++){
2284 for(j=0;j<3;j++) isen[j][i]=0;
2285 }
2286 nsen=0;
2287 for(i=0;i<3;i++){
2288 for(j=0;j<2;j++){
2289
2290 /***********************************************************
2291 (j=0) (j=1)
2292 (i=0) k=0,l=4,m=2+4=6 0->4->6 k=1,l=3,m=7 1->3->7
2293 (i=1) k=0,l=1,m=4+1=5 0->1->5 k=2,l=6,m=7 2->6->7
2294 (i=2) k=0,l=2,m=1+2=3 0->2->3 k=4,l=5,m=7 4->5->7
2295
2296 surface i-j direction
2297 0-0 -x
2298 0-1 +x
2299 1-0 -y
2300 1-1 +y
2301 2-0 -z
2302 2-1 +z
2303 ***********************************************************/
2304
2305 k=(int)(j*pow(2,i));
2306 l=(int)(k+pow(2,(i+2-j)%3));
2307 m=(int)(k+pow(2,(i+1)%3)+pow(2,(i+2)%3));
2308 x1=ten[l].x-ten[k].x;
2309 y1=ten[l].y-ten[k].y;
2310 z1=ten[l].z-ten[k].z;
2311 x2=ten[m].x-ten[l].x;
2312 y2=ten[m].y-ten[l].y;
2313 z2=ten[m].z-ten[l].z;
2314
2315 if((y1*z2-z1*y2)*(po->x-ten[l].x)+(z1*x2-x1*z2)*(po->y
2316 -ten[l].y)+(x1*y2-y1*x2)*(po->z-ten[l].z)>=0.){
2317
2318 /* if the plane is facing to the observer's eye do : */
2319 m=k; /* initialize present point */
2320 for(n1=0;n1<4;n1++){ // do for four sorrounding points
2321 l=m; /* define previous point */
2322
2323 /* redefine present point */
2324 if(n1==0) m=(int)(k+pow(2,(i+2-j)%3));
2325 if(n1==1) m=(int)(k+pow(2,(i+1)%3)+pow(2,(i+2)%3));
2326 if(n1==2) m=(int)(k+pow(2,(i+j+1)%3));
2327 if(n1==3) m=k;
2328 // now, the side start at point l and ends at point m
2329 if(nsen==0) brnch=1;
2330 else {
2331 brnch=1;
2332 for(n=0;n<nsen;n++){
2333 if(isen[1][n]==l && isen[0][n]==m){ // already there
2334 isen[2][n]=2; /* double line */
2335 brnch=0; // don't add this line
2336 break;
2337 }
2338 }
2339 }
2340 if(brnch==1){ // newly added line
2341 isen[0][nsen]=l; // start point
2342 isen[1][nsen]=m; // end point
2343 isen[2][nsen]=1; // single line
2344 nsen+=1;
2345 }
2346 } /* n1 */
2347 } /* facing this direction */
2348 } /* j */
2349 } /* i */
2350 blk->start=0; /* start line of a block */
2351 cont->start=0; /* start line of a contour */
2352
2353 /* first, collect single lines and place in the first part */
2354 line=0;
2355 for(i=0;i<nsen;i++){
2356 if(isen[2][i]==1){ // single line
2357 for(j=0;j<2;j++){
2358 k=isen[j][i];
2359 s.x=ten[k].x;
2360 s.y=ten[k].y;
2361 s.z=ten[k].z;
2362 FrameToBodyTAM(&s,po,&d);
2363 NormPoint3d(&d,&s,&r);
2364 ToPolar3d(&s,&direc);
2365 if(j==0){
2366 seg[line].sx=direc.x;
2367 seg[line].sy=1.5707963267-direc.y;
2368 }
2369 else{
2370 seg[line].ex=direc.x;
2371 seg[line].ey=1.5707963267-direc.y;
2372 }
2373 }
2374 line+=1;
2375 }
2376 }
2377
2378 cont->pte=line; /* next-to-end line of the contour */
2379
2380 /* then, collect double lines and place in the last part */
2381 for(i=0;i<nsen;i++){
2382 if(isen[2][i]==2){
2383 for(j=0;j<2;j++){
2384 k=isen[j][i];
2385 s.x=ten[k].x;
2386 s.y=ten[k].y;
2387 s.z=ten[k].z;
2388 FrameToBodyTAM(&s,po,&d);
2389 NormPoint3d(&d,&s,&r);
2390 ToPolar3d(&s,&direc);
2391 if(j==0){
2392 seg[line].sx=direc.x;
2393 seg[line].sy=1.5707963267-direc.y;
2394 }
2395 else{
2396 seg[line].ex=direc.x;
2397 seg[line].ey=1.5707963267-direc.y;
2398 }
2399 }
2400 line+=1;
2401 }
2402 }
2403 blk->pte=line; /* next-to-end line of the block */
2404}
2405
2406
2407/******************************************************************************
2408 FieldOfView
2409
2410 Get horizontal and vertical field of view to cover the object
2411******************************************************************************/
2412
2413void CMPbox::FieldOfView(TransAngleMatrix *vu,Point2d *az,Point2d *el)
2414{
2415 double x1,y1,z1,d;
2416 TransAngleEuler tae;
2417 TransAngleMatrix tam;
2418 int i;
2419 Point3d s;
2420 Point3d f;
2421 Point2d a;
2422
2423 x1=0.5*lx;
2424 y1=0.5*ly;
2425 z1=0.5*lz;
2426 tae.x=xg; // Translate Angle Euler for the part
2427 tae.y=yg;
2428 tae.z=zg;
2429 tae.theta=thetag;
2430 tae.phi=phig;
2431 tae.psi=psig;
2432 TAEtoTAM(&tae,&tam);
2433 az->x=10.;
2434 az->y=-10.;
2435 el->x=10.;
2436 el->y=-10.;
2437 for(i=0;i<8;i++){
2438 s.x=i&1? x1:-x1;
2439 s.y=i&2? y1:-y1;
2440 s.z=i&4? z1:-z1;
2441 BodyToFrameTAM(&s,&tam,&f);
2442 FrameToBodyTAM(&f,vu,&s);
2443 NormPoint3d(&s,&s,&d);
2444 ToPolar3d(&s,&a);
2445 if(a.x<az->x) az->x=a.x;
2446 if(a.x>az->y) az->y=a.x;
2447 if(a.y<el->x) el->x=a.y;
2448 if(a.y>el->y) el->y=a.y;
2449 }
2450}
2451
2452
2453/******************************************************************************
2454 CMPsphere - sphere
2455******************************************************************************/
2456
2457void CMPsphere::SetRadius(double x)
2458{
2459 r=x;
2460}
2461
2462
2463/******************************************************************************
2464 GetSurface
2465
2466 Return a line segment
2467
2468******************************************************************************/
2469
2470void CMPsphere::GetSurface(TransAngleMatrix *po,Surface *surf,Point2d *surfPt)
2471{
2472 Line3d eye;
2473 Point3d u;
2474 double d,ra,d1,r1;
2475 Point2d a;
2476 int i;
2477 double ang;
2478
2479 eye.sx=xg; // center of the sphere
2480 eye.sy=yg;
2481 eye.sz=zg;
2482 eye.ex=po->x; // thePolar of CThreeD - position of the viewer's eye
2483 eye.ey=po->y;
2484 eye.ez=po->z;
2485
2486 NormLine3d(&eye,&u,&d);
2487 ToPolar3d(&u,&a);
2488
2489 ra=r; // radius of the sphere
2490 d1=ra*ra/d; // distance to the center of tangent circle from sphere center
2491 r1=sqrt(ra*ra-d1*d1); // radius of the tangent circle
2492 r1*=1.03528; // r1/cos(15deg) to enable establishing the far-near
2493 // relationship with proximate parts.
2494
2495 surf->x0=xg+d1*u.x;
2496 surf->y0=yg+d1*u.y;
2497 surf->z0=zg+d1*u.z;
2498 surf->theta0=1.5707963268-a.y; // z axis toward viewer's eye
2499 surf->phi0=a.x+1.5707963268;
2500 surf->psi0=0.;
2501 surf->start=0;
2502 surf->pte=12;
2503 for(i=0;i<12;i++){
2504 ang=0.52359877559*i; // every 30 degree
2505 surfPt[i].x=r1*cos(ang);
2506 surfPt[i].y=r1*sin(ang);
2507 }
2508}
2509
2510
2511/******************************************************************************
2512 GetDrawData
2513
2514 Return data for drawing
2515
2516 po thePolar
2517
2518 seg array of line segments (clockwise angle from top, angle from line of sight)
2519 blk indices of segs which are to be drawn with lines
2520 cont indices of closed line segments. Inside of the closure
2521 will be painted by a color.
2522
2523******************************************************************************/
2524
2525void CMPsphere::GetDrawData(TransAngleMatrix *po,Line2d *seg,StartPTE *blk,StartPTE *cont)
2526{
2527 TransAngleEuler tae;
2528 TransAngleMatrix tam;
2529 Line3d eye;
2530 Point3d u,b,c;
2531 Point2d e,f,g;
2532 double d,ra,d1,r1;
2533 Point2d a;
2534 int i;
2535 double ang;
2536
2537 eye.sx=xg; // center of the sphere
2538 eye.sy=yg;
2539 eye.sz=zg;
2540 eye.ex=po->x; // thePolar of CThreeD - position of the viewer's eye
2541 eye.ey=po->y;
2542 eye.ez=po->z;
2543
2544 NormLine3d(&eye,&u,&d);
2545 ToPolar3d(&u,&a);
2546
2547 ra=r; // radius of the sphere
2548 d1=ra*ra/d; // distance to the center of tangent circle from sphere center
2549 r1=sqrt(ra*ra-d1*d1); // radius of the tangent circle
2550
2551 tae.x=xg+d1*u.x;
2552 tae.y=yg+d1*u.y;
2553 tae.z=zg+d1*u.z;
2554 tae.theta=1.5707963268-a.y; // z axis toward viewer's eye
2555 tae.phi=a.x+1.5707963268;
2556 tae.psi=0.;
2557 TAEtoTAM(&tae,&tam);
2558
2559 u.x=r1*cos(0.);
2560 u.y=r1*sin(0.);
2561 u.z=0.;
2562 BodyToFrameTAM(&u,&tam,&b);
2563 FrameToBodyTAM(&b,po,&u);
2564 NormPoint3d(&u,&b,&d);
2565 ToPolar3d(&b,&e);
2566 e.y=1.5707963267-e.y;
2567 f=e;
2568 for(i=0;i<90;i++){
2569 g=f;
2570 if(i<89){
2571 ang=.06981317008*(i+1); // every 4 degree
2572 u.x=r1*cos(ang);
2573 u.y=r1*sin(ang);
2574 u.z=0.;
2575 BodyToFrameTAM(&u,&tam,&b);
2576 FrameToBodyTAM(&b,po,&u);
2577 NormPoint3d(&u,&b,&d);
2578 ToPolar3d(&b,&f);
2579 f.y=1.5707963267-f.y;
2580 }
2581 else f=e;
2582 seg[i].sx=g.x;
2583 seg[i].sy=g.y;
2584 seg[i].ex=f.x;
2585 seg[i].ey=f.y;
2586
2587 }
2588 blk->start=0;
2589 blk->pte=90;
2590 cont->start=0;
2591 cont->pte=90;
2592}
2593
2594
2595
2596/******************************************************************************
2597 FieldOfView
2598
2599 Get horizontal and vertical field of view to cover the object
2600******************************************************************************/
2601
2602void CMPsphere::FieldOfView(TransAngleMatrix *vu,Point2d *az,Point2d *el)
2603{
2604 Point3d s;
2605 Point3d f;
2606 Point2d a;
2607 double b,d;
2608
2609 f.x=xg;
2610 f.y=yg;
2611 f.z=zg;
2612 FrameToBodyTAM(&f,vu,&s);
2613 NormPoint3d(&s,&s,&d);
2614 ToPolar3d(&s,&a);
2615 b=asin(r/d);
2616 az->x=a.x-b;
2617 az->y=a.x+b;
2618 el->x=a.y-b;
2619 el->y=a.y+b;
2620}
2621
2622/******************************************************************************
2623 CMPbodyOfRot - bodyOfRot
2624******************************************************************************/
2625
2626void CMPbodyOfRot::SetPoint(double x,double z)
2627{
2628 switch (npoint){
2629 case 0:
2630 p0.x=x;
2631 p0.y=z;
2632 npoint=1;
2633 break;
2634 case 1:
2635 p1.x=x;
2636 p1.y=z;
2637 npoint=2;
2638 break;
2639 default:
2640 break;
2641 }
2642}
2643
2644/******************************************************************************
2645 ClearVisSurface
2646
2647 Clear 'visibleSurface' instance variable so that new visibility
2648 from a new viewpoint can be obtained in GetSurface
2649******************************************************************************/
2650void CMPbodyOfRot::ClearVisSurface(void)
2651{
2652 visibleSurface=0;
2653}
2654
2655/******************************************************************************
2656 GetSurface
2657
2658 Return a surface
2659
2660 14 surfaces are related to an index number in the following way
2661
2662 index face toward...
2663 0 - 11 azimuth 30*index+15 (deg)
2664 12 +z
2665 13 -z
2666
2667 The instance variable visibleSurface is used to map between
2668 the above index and the reference index passed as the first parameter.
2669
2670 The returned bool value is false if no more visible surface exists
2671 for the asked reference number.
2672
2673******************************************************************************/
2674
2675bool CMPbodyOfRot::GetSurface(int n,bool nearside,TransAngleMatrix *po,
2676 Surface *surf,Point2d *surfPt)
2677{
2678 Point2d p2,p3;
2679 TransAngleEuler tae1;
2680 TransAngleMatrix tam1;
2681 double alpha,delta,r;
2682
2683 TransAngleEuler tae2;
2684 TransAngleEuler tae;
2685 Point3d s;
2686 Point3d e;
2687 int i,j,k,m;
2688
2689 p2.x=p0.x; // bottom radial point
2690 p2.y=p0.y;
2691 p3.x=p1.x; // top radial point
2692 p3.y=p1.y;
2693
2694 if(p2.x==0.) p2.x=p3.x*0.01; // cone
2695 if(p3.x==0.) p3.x=p2.x*0.01;
2696 if(p2.y==p3.y) { // didc
2697 if(p2.x>p3.x) p3.y+=0.01*p2.x;
2698 else p2.y+=0.01*p3.x;
2699 }
2700
2701 p2.x*=1.0178; // make the average radius same as the real radius
2702 p3.x*=1.0178;
2703
2704 tae1.x=xg; // Translate Angle Euler for the part
2705 tae1.y=yg;
2706 tae1.z=zg;
2707 tae1.theta=thetag;
2708 tae1.phi=phig;
2709 tae1.psi=psig;
2710
2711 if(visibleSurface==0){
2712 TAEtoTAM(&tae1,&tam1);
2713 j=1;
2714 delta=atan2((p2.x-p3.x)*cos(0.017453293*15.),p3.y-p2.y);
2715 r=-p2.y*(p3.x-p2.x)/(p3.y-p2.y)+p2.x;
2716 r*=cos(0.017453293*15.);
2717 s.z=r*cos(delta)*sin(delta);
2718 r=r*cos(delta)*cos(delta);
2719 for(i=0;i<12;i++){
2720 alpha=0.017453293*(30.*i+15.);
2721 s.x=r*cos(alpha);
2722 s.y=r*sin(alpha);
2723 BodyToFrameTAM(&s,&tam1,&e);
2724 if((xg-e.x)*(po->x-e.x)+(yg-e.y)*(po->y-e.y)
2725 +(zg-e.z)*(po->z-e.z)<0.) visibleSurface|=j;
2726 j=j<<1;
2727 }
2728 FrameToBodyTAM((Point3d *)po,&tam1,&s); // position of viewer
2729 if((s.z-p3.y)*(p2.y-p3.y)<0.) visibleSurface|=j; // upper surface
2730 j=j<<1;
2731 if((s.z-p2.y)*(p3.y-p2.y)<0.) visibleSurface|=j; // lower surface
2732
2733 // determine instance variable 'appear'
2734 j=1<<12;
2735 if(visibleSurface&j){
2736 j=1;
2737 k=0;
2738 for(i=0;i<12;i++){
2739 if(visibleSurface&j) k++;
2740 j=j<<1;
2741 }
2742 if(k==0) appear=1;
2743 else if(k==12) appear=3;
2744 else appear=2;
2745 }
2746 else{
2747 j=j<<1;
2748 if(visibleSurface&j){
2749 j=1;
2750 k=0;
2751 for(i=0;i<12;i++){
2752 if(visibleSurface&j) k++;
2753 j=j<<1;
2754 }
2755 if(k==0) appear=4;
2756 else if(k==12) appear=6;
2757 else appear=5;
2758 }
2759 else appear=0;
2760 }
2761 }
2762 if(isTube){
2763 if(appear==0||((appear==2||appear==5)&&nearside)){
2764 j=1;
2765 i=0;
2766 for(m=0;m<12;m++){
2767 if(visibleSurface&j){
2768 if(i==n) break;
2769 i++;
2770 }
2771 j=j<<1;
2772 }
2773 if(m==12) return false;
2774
2775 delta=atan2((p2.x-p3.x)*cos(0.017453293*15.),p3.y-p2.y);
2776 alpha=0.017453293*(30.*m+15.);
2777 r=p2.x*cos(0.017453293*15.); // radius to the center of the bottom segment
2778
2779 tae2.x=r*cos(alpha);
2780 tae2.y=r*sin(alpha);
2781 tae2.z=p2.y;
2782 tae2.theta=1.5707963268-delta;
2783 tae2.phi=alpha+1.5707963268;
2784 tae2.psi=0.;
2785
2786 surfPt[0].x=-r*tan(0.017453293*15.); // left-bottom
2787 surfPt[0].y=0.;
2788 surfPt[1].x=-surfPt[0].x; // right-bottom
2789 surfPt[1].y=0.;
2790 r=p3.x*cos(0.017453293*15.); // radius to the center of the top segment
2791 surfPt[2].x=r*tan(0.017453293*15.); // right-top
2792 surfPt[2].y=(p3.y-p2.y)/cos(delta);
2793 surfPt[3].x=-surfPt[2].x; // left-top
2794 surfPt[3].y=surfPt[2].y;
2795
2796 MoveTwiceTAE(&tae1,&tae2,&tae);
2797 surf->x0=tae.x;
2798 surf->y0=tae.y;
2799 surf->z0=tae.z;
2800 surf->theta0=tae.theta;
2801 surf->phi0=tae.phi;
2802 surf->psi0=tae.psi;
2803 surf->start=0;
2804 if(m<12) surf->pte=4;
2805 else surf->pte=12;
2806 return true;
2807 }
2808 if((appear==2||appear==5)&&!nearside){
2809 j=1;
2810 i=0;
2811 for(m=0;m<12;m++){
2812 if(!(visibleSurface&j)){
2813 if(i==n) break;
2814 i++;
2815 }
2816 j=j<<1;
2817 }
2818 if(m==12) return false;
2819
2820 delta=atan2((p2.x-p3.x)*cos(0.017453293*15.),p3.y-p2.y);
2821 alpha=0.017453293*(30.*m+15.);
2822 r=p2.x*cos(0.017453293*15.); // radius to the center of the bottom segment
2823
2824 tae2.x=r*cos(alpha);
2825 tae2.y=r*sin(alpha);
2826 tae2.z=p2.y;
2827 tae2.theta=1.5707963268-delta;
2828 tae2.phi=alpha+1.5707963268;
2829 tae2.psi=0.;
2830
2831 surfPt[0].x=r*tan(0.017453293*15.); // left-bottom
2832 surfPt[0].y=0.;
2833 surfPt[1].x=-surfPt[0].x; // right-bottom
2834 surfPt[1].y=0.;
2835 r=p3.x*cos(0.017453293*15.); // radius to the center of the top segment
2836 surfPt[2].x=-r*tan(0.017453293*15.); // right-top
2837 surfPt[2].y=(p3.y-p2.y)/cos(delta);
2838 surfPt[3].x=-surfPt[2].x; // left-top
2839 surfPt[3].y=surfPt[2].y;
2840
2841 MoveTwiceTAE(&tae1,&tae2,&tae);
2842 surf->x0=tae.x;
2843 surf->y0=tae.y;
2844 surf->z0=tae.z;
2845 surf->theta0=tae.theta;
2846 surf->phi0=tae.phi;
2847 surf->psi0=tae.psi;
2848 surf->start=0;
2849 if(m<12) surf->pte=4;
2850 else surf->pte=12;
2851 return true;
2852 }
2853 if(appear==1||appear==4||appear==3||appear==6){
2854 if(n>5) return false;
2855 if(nearside) m=n;
2856 else m=n+6;
2857
2858 delta=atan2((p2.x-p3.x)*cos(0.017453293*15.),p3.y-p2.y);
2859 alpha=0.017453293*(30.*m+15.);
2860 r=p2.x*cos(0.017453293*15.); // radius to the center of the bottom segment
2861 if(appear==1||appear==4){
2862 alpha-=3.14159265;
2863 delta=-delta;
2864 }
2865 tae2.x=r*cos(alpha);
2866 tae2.y=r*sin(alpha);
2867 tae2.z=p2.y;
2868 tae2.theta=1.5707963268-delta;
2869 tae2.phi=alpha+1.5707963268;
2870 tae2.psi=0.;
2871
2872 surfPt[0].x=-r*tan(0.017453293*15.); // left-bottom
2873 surfPt[0].y=0.;
2874 surfPt[1].x=-surfPt[0].x; // right-bottom
2875 surfPt[1].y=0.;
2876 r=p3.x*cos(0.017453293*15.); // radius to the center of the top segment
2877 surfPt[2].x=r*tan(0.017453293*15.); // right-top
2878 surfPt[2].y=(p3.y-p2.y)/cos(delta);
2879 surfPt[3].x=-surfPt[2].x; // left-top
2880 surfPt[3].y=surfPt[2].y;
2881
2882 MoveTwiceTAE(&tae1,&tae2,&tae);
2883 surf->x0=tae.x;
2884 surf->y0=tae.y;
2885 surf->z0=tae.z;
2886 surf->theta0=tae.theta;
2887 surf->phi0=tae.phi;
2888 surf->psi0=tae.psi;
2889 surf->start=0;
2890 if(m<12) surf->pte=4;
2891 else surf->pte=12;
2892 return true;
2893 }
2894 }
2895 else{ // solid
2896 j=1;
2897 i=0;
2898 for(m=0;m<14;m++){
2899 if(visibleSurface&j){
2900 if(i==n) break;
2901 i++;
2902 }
2903 j=j<<1;
2904 }
2905 if(m==14) return false;
2906
2907 if(m<12){
2908 delta=atan2((p2.x-p3.x)*cos(0.017453293*15.),p3.y-p2.y);
2909 alpha=0.017453293*(30.*m+15.);
2910 r=p2.x*cos(0.017453293*15.); // radius to the center of the bottom segment
2911
2912 tae2.x=r*cos(alpha);
2913 tae2.y=r*sin(alpha);
2914 tae2.z=p2.y;
2915 tae2.theta=1.5707963268-delta;
2916 tae2.phi=alpha+1.5707963268;
2917 tae2.psi=0.;
2918
2919 surfPt[0].x=-r*tan(0.017453293*15.); // left-bottom
2920 surfPt[0].y=0.;
2921 surfPt[1].x=-surfPt[0].x; // right-bottom
2922 surfPt[1].y=0.;
2923 r=p3.x*cos(0.017453293*15.); // radius to the center of the top segment
2924 surfPt[2].x=r*tan(0.017453293*15.); // right-top
2925 surfPt[2].y=(p3.y-p2.y)/cos(delta);
2926 surfPt[3].x=-surfPt[2].x; // left-top
2927 surfPt[3].y=surfPt[2].y;
2928 }
2929 if(m==12){
2930 tae2.x=0.;
2931 tae2.y=0.;
2932 tae2.z=p3.y;
2933 tae2.theta=0.;
2934 tae2.phi=0.;
2935 tae2.psi=0.;
2936
2937 r=p3.x;
2938 for(i=0;i<12;i++){
2939 alpha=0.017453293*30.*i;
2940 surfPt[i].x=r*cos(alpha);
2941 surfPt[i].y=r*sin(alpha);
2942 }
2943 }
2944 if(m==13){
2945 tae2.x=0.;
2946 tae2.y=0.;
2947 tae2.z=p2.y;
2948 tae2.theta=3.1415926536;
2949 tae2.phi=0.;
2950 tae2.psi=0.;
2951
2952 r=p2.x;
2953 for(i=0;i<12;i++){
2954 alpha=0.017453293*30.*i;
2955 surfPt[i].x=r*cos(alpha);
2956 surfPt[i].y=r*sin(alpha);
2957 }
2958 }
2959
2960 MoveTwiceTAE(&tae1,&tae2,&tae);
2961 surf->x0=tae.x;
2962 surf->y0=tae.y;
2963 surf->z0=tae.z;
2964 surf->theta0=tae.theta;
2965 surf->phi0=tae.phi;
2966 surf->psi0=tae.psi;
2967 surf->start=0;
2968 if(m<12) surf->pte=4;
2969 else surf->pte=12;
2970 return true;
2971 }
2972}
2973
2974/******************************************************************************
2975 GetCenter
2976
2977 Return the center point in frame coordinate
2978
2979******************************************************************************/
2980
2981void CMPbodyOfRot::GetCenter(Point3d *p)
2982{
2983 Point3d pt;
2984 TransAngleEuler tae;
2985 TransAngleMatrix tam;
2986
2987 pt.x=0.;
2988 pt.y=0.;
2989 pt.z=0.5*(p0.y+p1.y);
2990
2991 tae.x=xg;
2992 tae.y=yg;
2993 tae.z=zg;
2994 tae.theta=thetag;
2995 tae.phi=phig;
2996 tae.psi=psig;
2997
2998 TAEtoTAM(&tae,&tam);
2999 BodyToFrameTAM(&pt,&tam,p);
3000}
3001
3002/******************************************************************************
3003 GetDrawData
3004
3005 Return data for drawing
3006
3007 m segment number (0 - near side, 1 - far side)
3008 po thePolar
3009
3010 seg array of line segments (clockwise angle from top,
3011 angle from line of sight)
3012 blk indices of segs which are to be drawn with lines
3013 cont indices of closed line segments. Inside of the closure
3014 will be painted by a color.
3015
3016******************************************************************************/
3017
3018void CMPbodyOfRot::GetDrawData(int m,TransAngleMatrix *po,
3019 Line2d *seg,StartPTE *blk,StartPTE *cont)
3020{
3021 double eye[3];
3022 Point3d *ten0;
3023 Point3d *ten1;
3024 double x1,z1,x2,z2,y1,y2;
3025 TransAngleEuler tae;
3026 TransAngleMatrix tam1;
3027 int nseg,i,n0,n1,n2,j,k,is,is1,line,istart,iend,itb,ibt,n;
3028 double deld,a,co,si,r;
3029 Point3d q,s,d;
3030 Point2d direc;
3031
3032 if(oblong){
3033 GetDrawData1(po,seg,blk,cont);
3034 return;
3035 }
3036
3037 eye[0]=po->x;
3038 eye[1]=po->y;
3039 eye[2]=po->z;
3040
3041 // lower surface ten0, upper surface ten1
3042 // indices are from 0 to 90, 90th data being the same as 0th data
3043 ten0=new Point3d[91];
3044 ten1=new Point3d[91];
3045 x1=p0.x;
3046 z1=p0.y;
3047 x2=p1.x;
3048 z2=p1.y;
3049
3050 if(x1==0.) x1=x2*0.01; // cone
3051 if(x2==0.) x2=x1*0.01;
3052 if(z1==z2) { // didc
3053 if(x1>x2) z2+=0.01*x1;
3054 else z1-=0.01*x2;
3055 }
3056
3057 // set conversion parameters
3058 tae.x=xg;
3059 tae.y=yg;
3060 tae.z=zg;
3061 tae.theta=thetag;
3062 tae.phi=phig;
3063 tae.psi=psig;
3064
3065 TAEtoTAM(&tae,&tam1);
3066
3067 // At this point, initial conditions are set.
3068 // tam1 defines a local coordinate.
3069 // In the local coordinate, two points (x1,z1) and (x2,z2) are given.
3070 // The body of rotation is given by rotating
3071 // the line connecting the two points around z-axis.
3072 // The viewer is given by eye[]
3073
3074 /* fill coordinates of 2*nseg apexes */
3075
3076 nseg=90;
3077 deld=nseg;
3078 deld=2.*3.1415926536/deld;
3079 for(i=0;i<=nseg;i++){ // 91 data are set to each of the arrays
3080 a=i*deld;
3081 co=cos(a);
3082 si=sin(a);
3083 q.x=x1*co;
3084 q.y=x1*si;
3085 q.z=z1; /* lower surface */
3086 BodyToFrameTAM(&q,&tam1,&ten0[i]);
3087 q.x=x2*co;
3088 q.y=x2*si;
3089 q.z=z2; /* uppser surface */
3090 BodyToFrameTAM(&q,&tam1,&ten1[i]);
3091 }
3092
3093 // the outline was replaced by two circles ten0[] and ten1[]
3094 /* obtain side surface segments which are visible */
3095 n0=0;
3096 n1=0; // invisible-to-visible index
3097 n2=0; // visible-to-invisible index
3098
3099 // In CMPbodyOfRot::GetSurface, we have already categorized
3100 // into seven cases, which are inducated by the variable 'appear'.
3101 // Of those seven cases, the following three cases
3102 // 0 only the side is visible
3103 // 2 side and upper surface are visible
3104 // 5 side and lower surface are visible
3105 // have indices which separate visible and invisible portions of the side
3106 // surface.
3107
3108 if(appear==0||appear==2||appear==5){
3109 for(i=0;(i<=nseg)&&(n0!=2);i++){
3110
3111 /* each surface contains i-th and (i-1)th apexes */
3112 /* visibility of two adjacent surfaces are compared to */
3113 /* obtain both visible-to-invisible and invisible-to-visible */
3114 /* trasition indices */
3115
3116 j=i;
3117 if(j>=nseg) j-=nseg;
3118 k=j-1;
3119 if(k<0) k+=nseg;
3120 is=0;
3121 /* vector along lower surface */
3122 x1=ten0[j].x-ten0[k].x;
3123 y1=ten0[j].y-ten0[k].y;
3124 z1=ten0[j].z-ten0[k].z;
3125 /* vector directed from lower surface to upper sureace */
3126 x2=ten1[k].x-ten0[k].x;
3127 y2=ten1[k].y-ten0[k].y;
3128 z2=ten1[k].z-ten0[k].z;
3129 /* if surface is visible is = 1 */
3130 if((y1*z2-z1*y2)*(eye[0]-ten0[k].x)
3131 +(z1*x2-x1*z2)*(eye[1]-ten0[k].y)
3132 +(x1*y2-y1*x2)*(eye[2]-ten0[k].z)>0.) is=1; // visible
3133 if((i!=0)&&(is!=is1)) {
3134 if(is>is1) n1=k; /* k is invisible to visible */
3135 if(is<is1) n2=k; /* k is visible to invisible */
3136 n0++;
3137 }
3138 if(n0!=2) is1=is;
3139 }
3140 }
3141
3142 /* now, categorization is over, process each cases */
3143
3144 // The order of processing is as followes
3145 // appear solid tube-near tube-far
3146 // 0 x x x
3147 // 2,5 x
3148 // 2,5 x
3149 // 5 x
3150 // 2 x
3151 // 1,3,4,6 x
3152 // 1,3,4,6 x
3153 // 1,3,4,6 x
3154
3155 line=0;
3156 if(appear==0||((appear==2||appear==5)&&isTube&&m==0)){
3157 // appear=0,both for a solid and a tube, near and far side
3158 // appear=2 or 5, near side of the tube
3159 /* side only */
3160 blk->start=line;
3161 cont->start=line;
3162 n0=n2-n1; /* number of visible points - 1 */
3163 if(n0<0) n0+=nseg;
3164 for(i=0;i<n0;i++){
3165 j=n1+i+1;
3166 if(j>=nseg) j-=nseg;
3167 k=j-1;
3168 if(k<0) k+=nseg;
3169 /* along lower surface normal direction (to the order of smaller to larger indices) */
3170 s=ten0[k];
3171 FrameToBodyTAM(&s,po,&d);
3172 NormPoint3d(&d,&s,&r);
3173 ToPolar3d(&s,&direc);
3174 seg[line].sx=direc.x;
3175 seg[line].sy=1.5707963267-direc.y;
3176 s=ten0[j];
3177 FrameToBodyTAM(&s,po,&d);
3178 NormPoint3d(&d,&s,&r);
3179 ToPolar3d(&s,&direc);
3180 seg[line].ex=direc.x;
3181 seg[line].ey=1.5707963267-direc.y;
3182 line++;
3183 /* along upper surface reversed order */
3184 s=ten1[j];
3185 FrameToBodyTAM(&s,po,&d);
3186 NormPoint3d(&d,&s,&r);
3187 ToPolar3d(&s,&direc);
3188 seg[line].sx=direc.x;
3189 seg[line].sy=1.5707963267-direc.y;
3190 s=ten1[k];
3191 FrameToBodyTAM(&s,po,&d);
3192 NormPoint3d(&d,&s,&r);
3193 ToPolar3d(&s,&direc);
3194 seg[line].ex=direc.x;
3195 seg[line].ey=1.5707963267-direc.y;
3196 line++;
3197 }
3198 /* invisible-to-visible side line from top to bottom */
3199 s=ten1[n1];
3200 FrameToBodyTAM(&s,po,&d);
3201 NormPoint3d(&d,&s,&r);
3202 ToPolar3d(&s,&direc);
3203 seg[line].sx=direc.x;
3204 seg[line].sy=1.5707963267-direc.y;
3205 s=ten0[n1];
3206 FrameToBodyTAM(&s,po,&d);
3207 NormPoint3d(&d,&s,&r);
3208 ToPolar3d(&s,&direc);
3209 seg[line].ex=direc.x;
3210 seg[line].ey=1.5707963267-direc.y;
3211 line++;
3212 /* visible-to-invisible side line from bottom to top */
3213 s=ten0[n2];
3214 FrameToBodyTAM(&s,po,&d);
3215 NormPoint3d(&d,&s,&r);
3216 ToPolar3d(&s,&direc);
3217 seg[line].sx=direc.x;
3218 seg[line].sy=1.5707963267-direc.y;
3219 s=ten1[n2];
3220 FrameToBodyTAM(&s,po,&d);
3221 NormPoint3d(&d,&s,&r);
3222 ToPolar3d(&s,&direc);
3223 seg[line].ex=direc.x;
3224 seg[line].ey=1.5707963267-direc.y;
3225 line++;
3226 blk->pte=line;
3227 cont->pte=line;
3228 delete[] ten0; // DisposPtr((Ptr)ten0);
3229 delete[] ten1; // DisposPtr((Ptr)ten1);
3230 return;
3231 }
3232
3233 if((appear==2||appear==5)&&isTube&&m==1){
3234 // appear=2 or 5, far side of the tube
3235 /* side only */
3236 blk->start=line;
3237 cont->start=line;
3238 n0=n1-n2; /* number of invisible points - 1 */
3239 if(n0<0) n0+=nseg;
3240 for(i=0;i<n0;i++){
3241 j=n2+i+1;
3242 if(j>=nseg) j-=nseg;
3243 k=j-1;
3244 if(k<0) k+=nseg;
3245 /* along lower surface in reverse direction */
3246 s=ten0[j];
3247 FrameToBodyTAM(&s,po,&d);
3248 NormPoint3d(&d,&s,&r);
3249 ToPolar3d(&s,&direc);
3250 seg[line].sx=direc.x;
3251 seg[line].sy=1.5707963267-direc.y;
3252 s=ten0[k];
3253 FrameToBodyTAM(&s,po,&d);
3254 NormPoint3d(&d,&s,&r);
3255 ToPolar3d(&s,&direc);
3256 seg[line].ex=direc.x;
3257 seg[line].ey=1.5707963267-direc.y;
3258 line++;
3259 /* along upper surface in normal direction */
3260 s=ten1[k];
3261 FrameToBodyTAM(&s,po,&d);
3262 NormPoint3d(&d,&s,&r);
3263 ToPolar3d(&s,&direc);
3264 seg[line].sx=direc.x;
3265 seg[line].sy=1.5707963267-direc.y;
3266 s=ten1[j];
3267 FrameToBodyTAM(&s,po,&d);
3268 NormPoint3d(&d,&s,&r);
3269 ToPolar3d(&s,&direc);
3270 seg[line].ex=direc.x;
3271 seg[line].ey=1.5707963267-direc.y;
3272 line++;
3273 }
3274 /* invisible-to-visible side line from top tp bottom */
3275 s=ten1[n1];
3276 FrameToBodyTAM(&s,po,&d);
3277 NormPoint3d(&d,&s,&r);
3278 ToPolar3d(&s,&direc);
3279 seg[line].sx=direc.x;
3280 seg[line].sy=1.5707963267-direc.y;
3281 s=ten0[n1];
3282 FrameToBodyTAM(&s,po,&d);
3283 NormPoint3d(&d,&s,&r);
3284 ToPolar3d(&s,&direc);
3285 seg[line].ex=direc.x;
3286 seg[line].ey=1.5707963267-direc.y;
3287 line++;
3288 /* visible-to-invisible side line from bottom to top */
3289 s=ten0[n2];
3290 FrameToBodyTAM(&s,po,&d);
3291 NormPoint3d(&d,&s,&r);
3292 ToPolar3d(&s,&direc);
3293 seg[line].sx=direc.x;
3294 seg[line].sy=1.5707963267-direc.y;
3295 s=ten1[n2];
3296 FrameToBodyTAM(&s,po,&d);
3297 NormPoint3d(&d,&s,&r);
3298 ToPolar3d(&s,&direc);
3299 seg[line].ex=direc.x;
3300 seg[line].ey=1.5707963267-direc.y;
3301 line++;
3302 blk->pte=line;
3303 cont->pte=line;
3304 delete[] ten0; // DisposPtr((Ptr)ten0);
3305 delete[] ten1; // DisposPtr((Ptr)ten1);
3306 return;
3307 }
3308
3309 if(appear==5&&!isTube){ // for a solid
3310 /* side and lower surface */
3311 blk->start=line;
3312 cont->start=line;
3313 n0=n2-n1;
3314 if(n0<0) n0+=nseg;
3315 for(i=0;i<n0;i++){
3316 j=n2-i-1;
3317 if(j<0) j+=nseg;
3318 k=j+1;
3319 if(k>=nseg) k-=nseg;
3320 /* contor along the upper surface */
3321 s=ten1[k];
3322 FrameToBodyTAM(&s,po,&d);
3323 NormPoint3d(&d,&s,&r);
3324 ToPolar3d(&s,&direc);
3325 seg[line].sx=direc.x;
3326 seg[line].sy=1.5707963267-direc.y;
3327 s=ten1[j];
3328 FrameToBodyTAM(&s,po,&d);
3329 NormPoint3d(&d,&s,&r);
3330 ToPolar3d(&s,&direc);
3331 seg[line].ex=direc.x;
3332 seg[line].ey=1.5707963267-direc.y;
3333 line++;
3334 }
3335 s=ten1[n1];
3336 FrameToBodyTAM(&s,po,&d);
3337 NormPoint3d(&d,&s,&r);
3338 ToPolar3d(&s,&direc);
3339 seg[line].sx=direc.x;
3340 seg[line].sy=1.5707963267-direc.y;
3341 s=ten0[n1];
3342 FrameToBodyTAM(&s,po,&d);
3343 NormPoint3d(&d,&s,&r);
3344 ToPolar3d(&s,&direc);
3345 seg[line].ex=direc.x;
3346 seg[line].ey=1.5707963267-direc.y;
3347 line++;
3348 s=ten0[n2];
3349 FrameToBodyTAM(&s,po,&d);
3350 NormPoint3d(&d,&s,&r);
3351 ToPolar3d(&s,&direc);
3352 seg[line].sx=direc.x;
3353 seg[line].sy=1.5707963267-direc.y;
3354 s=ten1[n2];
3355 FrameToBodyTAM(&s,po,&d);
3356 NormPoint3d(&d,&s,&r);
3357 ToPolar3d(&s,&direc);
3358 seg[line].ex=direc.x;
3359 seg[line].ey=1.5707963267-direc.y;
3360 line++;
3361 n0=n1-n2;
3362 if(n0<0) n0=n0+nseg;
3363 cont->pte=line+n0;
3364 /* part of the followings constitute contour */
3365 for(i=0;i<nseg;i++){
3366 j=n2+i+1;
3367 if(j>=nseg) j=j-nseg;
3368 k=j-1;
3369 if(k<0) k+=nseg;
3370 s=ten0[j];
3371 FrameToBodyTAM(&s,po,&d);
3372 NormPoint3d(&d,&s,&r);
3373 ToPolar3d(&s,&direc);
3374 seg[line].sx=direc.x;
3375 seg[line].sy=1.5707963267-direc.y;
3376 s=ten0[k];
3377 FrameToBodyTAM(&s,po,&d);
3378 NormPoint3d(&d,&s,&r);
3379 ToPolar3d(&s,&direc);
3380 seg[line].ex=direc.x;
3381 seg[line].ey=1.5707963267-direc.y;
3382 line++;
3383 }
3384 blk->pte=line;
3385 delete[] ten0; // DisposPtr((Ptr)ten0);
3386 delete[] ten1; // DisposPtr((Ptr)ten1);
3387 return;
3388 }
3389 if(appear==2&&!isTube){ // for a solid
3390 /* side and upper surface */
3391 blk->start=line;
3392 cont->start=line;
3393 n0=n2-n1;
3394 if(n0<0) n0+=nseg;
3395 /* line along lower surface */
3396 for(i=0;i<n0;i++){
3397 j=n1+i+1;
3398 if(j>=nseg) j-=nseg;
3399 k=j-1;
3400 if(k<0) k+=nseg;
3401 s=ten0[k];
3402 FrameToBodyTAM(&s,po,&d);
3403 NormPoint3d(&d,&s,&r);
3404 ToPolar3d(&s,&direc);
3405 seg[line].sx=direc.x;
3406 seg[line].sy=1.5707963267-direc.y;
3407 s=ten0[j];
3408 FrameToBodyTAM(&s,po,&d);
3409 NormPoint3d(&d,&s,&r);
3410 ToPolar3d(&s,&direc);
3411 seg[line].ex=direc.x;
3412 seg[line].ey=1.5707963267-direc.y;
3413 line++;
3414 }
3415 s=ten1[n1];
3416 FrameToBodyTAM(&s,po,&d);
3417 NormPoint3d(&d,&s,&r);
3418 ToPolar3d(&s,&direc);
3419 seg[line].sx=direc.x;
3420 seg[line].sy=1.5707963267-direc.y;
3421 s=ten0[n1];
3422 FrameToBodyTAM(&s,po,&d);
3423 NormPoint3d(&d,&s,&r);
3424 ToPolar3d(&s,&direc);
3425 seg[line].ex=direc.x;
3426 seg[line].ey=1.5707963267-direc.y;
3427 line++;
3428 s=ten0[n2];
3429 FrameToBodyTAM(&s,po,&d);
3430 NormPoint3d(&d,&s,&r);
3431 ToPolar3d(&s,&direc);
3432 seg[line].sx=direc.x;
3433 seg[line].sy=1.5707963267-direc.y;
3434 s=ten1[n2];
3435 FrameToBodyTAM(&s,po,&d);
3436 NormPoint3d(&d,&s,&r);
3437 ToPolar3d(&s,&direc);
3438 seg[line].ex=direc.x;
3439 seg[line].ey=1.5707963267-direc.y;
3440 line++;
3441 n0=n1-n2;
3442 if(n0<0) n0=n0+nseg;
3443 cont->pte=line+n0;
3444 for(i=0;i<nseg;i++){
3445 j=n2+i+1;
3446 if(j>=nseg) j=j-nseg;
3447 k=j-1;
3448 if(k<0) k+=nseg;
3449 s=ten1[k];
3450 FrameToBodyTAM(&s,po,&d);
3451 NormPoint3d(&d,&s,&r);
3452 ToPolar3d(&s,&direc);
3453 seg[line].sx=direc.x;
3454 seg[line].sy=1.5707963267-direc.y;
3455 s=ten1[j];
3456 FrameToBodyTAM(&s,po,&d);
3457 NormPoint3d(&d,&s,&r);
3458 ToPolar3d(&s,&direc);
3459 seg[line].ex=direc.x;
3460 seg[line].ey=1.5707963267-direc.y;
3461 line++;
3462 }
3463 blk->pte=line;
3464 delete[] ten0; // DisposPtr((Ptr)ten0);
3465 delete[] ten1; // DisposPtr((Ptr)ten1);
3466 return;
3467 }
3468 if(!isTube){ // for a solid, appear=1,3,4 or 6
3469 /* lower or upper surface only */
3470 blk->start=line;
3471 cont->pte=line;
3472 for(i=0;i<nseg;i++){ // contour circle
3473 j=i;
3474 k=j-1;
3475 if(k<0) k+=nseg;
3476 if(appear==1||appear==6){ // contour is upper surface
3477 if(appear==1){
3478 n=k;
3479 k=j;
3480 j=n;
3481 }
3482 s=ten1[j];
3483 FrameToBodyTAM(&s,po,&d);
3484 NormPoint3d(&d,&s,&r);
3485 ToPolar3d(&s,&direc);
3486 seg[line].sx=direc.x;
3487 seg[line].sy=1.5707963267-direc.y;
3488 s=ten1[k];
3489 FrameToBodyTAM(&s,po,&d);
3490 NormPoint3d(&d,&s,&r);
3491 ToPolar3d(&s,&direc);
3492 seg[line].ex=direc.x;
3493 seg[line].ey=1.5707963267-direc.y;
3494 }
3495 else{ // appear=3 or 4
3496 if(appear==3){
3497 n=k;
3498 k=j;
3499 j=n;
3500 }
3501 s=ten0[j];
3502 FrameToBodyTAM(&s,po,&d);
3503 NormPoint3d(&d,&s,&r);
3504 ToPolar3d(&s,&direc);
3505 seg[line].sx=direc.x;
3506 seg[line].sy=1.5707963267-direc.y;
3507 s=ten0[k];
3508 FrameToBodyTAM(&s,po,&d);
3509 NormPoint3d(&d,&s,&r);
3510 ToPolar3d(&s,&direc);
3511 seg[line].ex=direc.x;
3512 seg[line].ey=1.5707963267-direc.y;
3513 }
3514 line++;
3515 }
3516 cont->pte=line;
3517 if(appear==3){
3518 for(i=0;i<nseg;i++){ // upper surface is inner circle
3519 j=i;
3520 k=j-1;
3521 if(k<0) k+=nseg;
3522 s=ten1[j];
3523 FrameToBodyTAM(&s,po,&d);
3524 NormPoint3d(&d,&s,&r);
3525 ToPolar3d(&s,&direc);
3526 seg[line].sx=direc.x;
3527 seg[line].sy=1.5707963267-direc.y;
3528 s=ten1[k];
3529 FrameToBodyTAM(&s,po,&d);
3530 NormPoint3d(&d,&s,&r);
3531 ToPolar3d(&s,&direc);
3532 seg[line].ex=direc.x;
3533 seg[line].ey=1.5707963267-direc.y;
3534 }
3535 line++;
3536 }
3537 if(appear==6){
3538 for(i=0;i<nseg;i++){ // lower surface is inner circle
3539 j=i;
3540 k=j-1;
3541 if(k<0) k+=nseg;
3542 s=ten0[j];
3543 FrameToBodyTAM(&s,po,&d);
3544 NormPoint3d(&d,&s,&r);
3545 ToPolar3d(&s,&direc);
3546 seg[line].sx=direc.x;
3547 seg[line].sy=1.5707963267-direc.y;
3548 s=ten0[k];
3549 FrameToBodyTAM(&s,po,&d);
3550 NormPoint3d(&d,&s,&r);
3551 ToPolar3d(&s,&direc);
3552 seg[line].ex=direc.x;
3553 seg[line].ey=1.5707963267-direc.y;
3554 }
3555 line++;
3556 }
3557 blk->pte=line;
3558 delete[] ten0; // DisposPtr((Ptr)ten0);
3559 delete[] ten1; // DisposPtr((Ptr)ten1);
3560 return;
3561 }
3562 // remaining cases are appear 1,3,4 and 6 of tubes
3563
3564 // The upper surface is given by an array ten1[].
3565 // The lower surface is given by an array ten0[].
3566 // looked from... outer circle inner circle
3567 // appear=1 the upper side upper circle (normal) lower circle (reverse)
3568 // appear=3 the upper side lower circle (normal) upper circle (reverse)
3569 // appear=4 the lower side lower circle (reverse) upper circle (normal)
3570 // appear=6 the lower side upper circle (reverse) lower circle (normal)
3571
3572 cont->start=line;
3573
3574 if(m==0){ // near side
3575 istart=0;
3576 iend=45;
3577 }
3578 else{ // far side
3579 istart=45;
3580 iend=90;
3581 }
3582 // 2 segments, which are part of the contour,
3583 // but, not part of the block.
3584
3585 if(appear==1||appear==4){
3586 itb=iend; // top to botton transition
3587 ibt=istart; // bottom to top transition
3588 }
3589 if(appear==3||appear==6){
3590 itb=istart; // top to bottom transition
3591 ibt=iend;
3592 }
3593 s=ten1[itb];
3594 FrameToBodyTAM(&s,po,&d);
3595 NormPoint3d(&d,&s,&r);
3596 ToPolar3d(&s,&direc);
3597 seg[line].sx=direc.x;
3598 seg[line].sy=1.5707963267-direc.y;
3599 s=ten0[itb];
3600 FrameToBodyTAM(&s,po,&d);
3601 NormPoint3d(&d,&s,&r);
3602 ToPolar3d(&s,&direc);
3603 seg[line].ex=direc.x;
3604 seg[line].ey=1.5707963267-direc.y;
3605 line++;
3606 s=ten0[ibt];
3607 FrameToBodyTAM(&s,po,&d);
3608 NormPoint3d(&d,&s,&r);
3609 ToPolar3d(&s,&direc);
3610 seg[line].sx=direc.x;
3611 seg[line].sy=1.5707963267-direc.y;
3612 s=ten1[ibt];
3613 FrameToBodyTAM(&s,po,&d);
3614 NormPoint3d(&d,&s,&r);
3615 ToPolar3d(&s,&direc);
3616 seg[line].ex=direc.x;
3617 seg[line].ey=1.5707963267-direc.y;
3618 line++;
3619
3620
3621 // now, the circular portion
3622 blk->start=line;
3623
3624 for(i=istart;i<iend;i++){
3625 if(appear==1||appear==4){
3626 j=i;
3627 k=i+1;
3628 }
3629 else{ // appear==3||appear==6
3630 j=i+1;
3631 k=i;
3632 }
3633 s=ten1[j];
3634 FrameToBodyTAM(&s,po,&d);
3635 NormPoint3d(&d,&s,&r);
3636 ToPolar3d(&s,&direc);
3637 seg[line].sx=direc.x;
3638 seg[line].sy=1.5707963267-direc.y;
3639 s=ten1[k];
3640 FrameToBodyTAM(&s,po,&d);
3641 NormPoint3d(&d,&s,&r);
3642 ToPolar3d(&s,&direc);
3643 seg[line].ex=direc.x;
3644 seg[line].ey=1.5707963267-direc.y;
3645 line++;
3646 }
3647 for(i=istart;i<iend;i++){
3648 if(appear==1||appear==4){
3649 j=i+1;
3650 k=i;
3651 }
3652 else{ // appear==3||appear==6
3653 j=i;
3654 k=i+1;
3655 }
3656 s=ten0[j];
3657 FrameToBodyTAM(&s,po,&d);
3658 NormPoint3d(&d,&s,&r);
3659 ToPolar3d(&s,&direc);
3660 seg[line].sx=direc.x;
3661 seg[line].sy=1.5707963267-direc.y;
3662 s=ten0[k];
3663 FrameToBodyTAM(&s,po,&d);
3664 NormPoint3d(&d,&s,&r);
3665 ToPolar3d(&s,&direc);
3666 seg[line].ex=direc.x;
3667 seg[line].ey=1.5707963267-direc.y;
3668 line++;
3669 }
3670 cont->pte=line;
3671 blk->pte=line;
3672 delete[] ten0; // DisposPtr((Ptr)ten0);
3673 delete[] ten1; // DisposPtr((Ptr)ten1);
3674}
3675
3676/******************************************************************************
3677 GetDrawData1
3678
3679 Return data for drawing for the case of oblong
3680
3681 po thePolar
3682
3683 seg array of line segments (clockwise angle from top,
3684 angle from line of sight)
3685 blk indices of segs which are to be drawn with lines
3686 cont indices of closed line segments. Inside of the closure
3687 will be painted by a color.
3688
3689******************************************************************************/
3690
3691void CMPbodyOfRot::GetDrawData1(TransAngleMatrix *po,
3692 Line2d *seg,StartPTE *blk,StartPTE *cont)
3693{
3694 int i,n;
3695 double z1,z2,r,ang0,ang1,ang2,ang,theta,dthet,d,d1,r1,d2,r2,l1,l2;
3696 TransAngleEuler tae;
3697 TransAngleMatrix tam;
3698 Point3d q,s1,u1,s2,u2,v,w,u,b;
3699 Line3d eye,tz,ty;
3700 Point2d e,f,g,a;
3701
3702 r=p0.x; // radius of sphere
3703 z1=p0.y+r; // bottom center of sphere
3704 z2=p1.y-r; // top center of sphere
3705
3706 tae.x=xg;
3707 tae.y=yg;
3708 tae.z=zg;
3709 tae.theta=thetag;
3710 tae.phi=phig;
3711 tae.psi=psig;
3712 TAEtoTAM(&tae,&tam);
3713
3714 q.x=0.;
3715 q.y=0.;
3716 q.z=z1; /* bottom center of sphere */
3717 BodyToFrameTAM(&q,&tam,&s1);
3718 eye.sx=s1.x; // center of the sphere
3719 eye.sy=s1.y;
3720 eye.sz=s1.z;
3721 eye.ex=po->x; // thePolar of CThreeD - position of the viewer's eye
3722 eye.ey=po->y;
3723 eye.ez=po->z;
3724
3725 NormLine3d(&eye,&u1,&d);
3726 d1=r*r/d; // distance to the center of tangent circle from sphere center
3727 r1=sqrt(r*r-d1*d1); // radius of the tangent circle
3728 l1=d-d1; // dist from circle center to eye
3729
3730 q.x=0.;
3731 q.y=0.;
3732 q.z=z2; /* top center of sphere */
3733 BodyToFrameTAM(&q,&tam,&s2);
3734 eye.sx=s2.x; // center of the sphere
3735 eye.sy=s2.y;
3736 eye.sz=s2.z;
3737 eye.ex=po->x; // thePolar of CThreeD - position of the viewer's eye
3738 eye.ey=po->y;
3739 eye.ez=po->z;
3740
3741 NormLine3d(&eye,&u2,&d);
3742 d2=r*r/d; // distance to the center of tangent circle from sphere center
3743 r2=sqrt(r*r-d2*d2); // radius of the tangent circle
3744 l2=d-d2; // dist from circle center to eye
3745
3746 ang0=acos(u1.x*u2.x+u1.y*u2.y+u1.z*u2.z); // angular separation of 2 centers
3747 ang1=atan2(r1,l1); // half-cone angle of bottom sphere
3748 ang2=atan2(r2,l2); // half-cone angle of top sphere
3749
3750 if(ang1>=0.98*(ang0+ang2)||ang2>=0.98*(ang0+ang1)){ // only one sphere is visible
3751 if(ang2>ang1){ // only top sphere is visible
3752 s1=s2; // sphere center position
3753 d1=d2; // dstance to the center of circle from sphere center
3754 r1=r2; // radius of tangent circle
3755 u1=u2; // unit vector toward viewer
3756 }
3757 ToPolar3d(&u1,&a);
3758
3759 tae.x=s1.x+d1*u1.x;
3760 tae.y=s1.y+d1*u1.y;
3761 tae.z=s1.z+d1*u1.z;
3762 tae.theta=1.5707963268-a.y; // z axis toward viewer's eye
3763 tae.phi=a.x+1.5707963268;
3764 tae.psi=0.;
3765 TAEtoTAM(&tae,&tam);
3766
3767 u.x=r1*cos(0.);
3768 u.y=r1*sin(0.);
3769 u.z=0.;
3770 BodyToFrameTAM(&u,&tam,&b);
3771 FrameToBodyTAM(&b,po,&u);
3772 NormPoint3d(&u,&b,&d);
3773 ToPolar3d(&b,&e);
3774 e.y=1.5707963267-e.y;
3775 f=e;
3776 for(i=0;i<90;i++){
3777 g=f;
3778 if(i<89){
3779 ang=.06981317008*(i+1); // every 4 degree
3780 u.x=r1*cos(ang);
3781 u.y=r1*sin(ang);
3782 u.z=0.;
3783 BodyToFrameTAM(&u,&tam,&b);
3784 FrameToBodyTAM(&b,po,&u);
3785 NormPoint3d(&u,&b,&d);
3786 ToPolar3d(&b,&f);
3787 f.y=1.5707963267-f.y;
3788 }
3789 else f=e;
3790 seg[i].sx=g.x;
3791 seg[i].sy=g.y;
3792 seg[i].ex=f.x;
3793 seg[i].ey=f.y;
3794
3795 }
3796 blk->start=0;
3797 blk->pte=90;
3798 cont->start=0;
3799 cont->pte=90;
3800 return;
3801 }
3802 // cylinder part is visible
3803 theta=1.5707963268+asin((ang1-ang2)/ang0); // half-visible angle of bottom circle
3804 // construct coordinate system with z toward eye and x toward bottom from top
3805 v.x=u2.x-u1.x; // unnormalized vector from top to bottom
3806 v.y=u2.y-u1.y;
3807 v.z=u2.z-u1.z;
3808 tz.sx=s1.x+d1*u1.x; // center of bottom circle
3809 tz.sy=s1.y+d1*u1.y;
3810 tz.sz=s1.z+d1*u1.z;
3811 tz.ex=tz.sx+u1.x; // vector toward viewer
3812 tz.ey=tz.sy+u1.y;
3813 tz.ez=tz.sz+u1.z;
3814
3815 VectorProduct(&u1,&v,&w); // obtain y-axis
3816 ty=tz;
3817 ty.ex=ty.sx+w.x; // vector toward viewer
3818 ty.ey=ty.sy+w.y;
3819 ty.ez=ty.sz+w.z;
3820 LinesToTAM(&ty,&tz,yzAxis,&tam);
3821
3822 dthet=.06981317008; // angle of a segment
3823 n=(int)(theta/dthet); // half number of segments
3824 theta=n*dthet;
3825 u.x=r1*cos(-theta);
3826 u.y=r1*sin(-theta);
3827 u.z=0.;
3828 BodyToFrameTAM(&u,&tam,&b);
3829 FrameToBodyTAM(&b,po,&u);
3830 NormPoint3d(&u,&b,&d);
3831 ToPolar3d(&b,&e);
3832 e.y=1.5707963267-e.y;
3833 f=e;
3834 n*=2;
3835 for(i=0;i<n;i++){
3836 g=f;
3837 ang=dthet*(i+1)-theta; // -theta+dthet to theta
3838 u.x=r1*cos(ang);
3839 u.y=r1*sin(ang);
3840 u.z=0.;
3841 BodyToFrameTAM(&u,&tam,&b);
3842 FrameToBodyTAM(&b,po,&u);
3843 NormPoint3d(&u,&b,&d);
3844 ToPolar3d(&b,&f);
3845 f.y=1.5707963267-f.y;
3846 seg[i].sx=g.x;
3847 seg[i].sy=g.y;
3848 seg[i].ex=f.x;
3849 seg[i].ey=f.y;
3850
3851 }
3852
3853 // construct coordinate system with z toward eye and x toward bottom from top
3854 v.x=u2.x-u1.x; // unnormalized vector from top to bottom
3855 v.y=u2.y-u1.y;
3856 v.z=u2.z-u1.z;
3857 tz.sx=s2.x+d2*u2.x; // center of top circle
3858 tz.sy=s2.y+d2*u2.y;
3859 tz.sz=s2.z+d2*u2.z;
3860 tz.ex=tz.sx+u2.x; // vector toward viewer
3861 tz.ey=tz.sy+u2.y;
3862 tz.ez=tz.sz+u2.z;
3863
3864 VectorProduct(&u2,&v,&w); // obtain y-axis
3865 ty=tz;
3866 ty.ex=ty.sx+w.x; // vector toward viewer
3867 ty.ey=ty.sy+w.y;
3868 ty.ez=ty.sz+w.z;
3869 LinesToTAM(&ty,&tz,yzAxis,&tam);
3870
3871 for(;i<91;i++){
3872 g=f;
3873 ang=dthet*i-theta;
3874 u.x=r2*cos(ang);
3875 u.y=r2*sin(ang);
3876 u.z=0.;
3877 BodyToFrameTAM(&u,&tam,&b);
3878 FrameToBodyTAM(&b,po,&u);
3879 NormPoint3d(&u,&b,&d);
3880 ToPolar3d(&b,&f);
3881 f.y=1.5707963267-f.y;
3882 seg[i].sx=g.x;
3883 seg[i].sy=g.y;
3884 seg[i].ex=f.x;
3885 seg[i].ey=f.y;
3886 }
3887 seg[i].sx=f.x;
3888 seg[i].sy=f.y;
3889 seg[i].ex=e.x;
3890 seg[i].ey=e.y;
3891 blk->start=0;
3892 blk->pte=92;
3893 cont->start=0;
3894 cont->pte=92;
3895}
3896
3897
3898/******************************************************************************
3899 FieldOfView
3900
3901 Get horizontal and vertical field of view to cover the object
3902******************************************************************************/
3903
3904void CMPbodyOfRot::FieldOfView(TransAngleMatrix *vu,Point2d *az,Point2d *el)
3905{
3906 TransAngleEuler tae;
3907 TransAngleMatrix tam;
3908 int i;
3909 Point3d s;
3910 Point3d f;
3911 Point2d a;
3912 double r0,z0,r1,z1,ca,sa,ang,d;
3913
3914 r0=p0.x;
3915 z0=p0.y;
3916 r1=p1.x;
3917 z1=p1.y;
3918 tae.x=xg; // Translate Angle Euler for the part
3919 tae.y=yg;
3920 tae.z=zg;
3921 tae.theta=thetag;
3922 tae.phi=phig;
3923 tae.psi=psig;
3924 TAEtoTAM(&tae,&tam);
3925 az->x=10.;
3926 az->y=-10.;
3927 el->x=10.;
3928 el->y=-10.;
3929 for(i=0;i<12;i++){
3930 ang=0.01745329*30.*i;
3931 ca=cos(ang);
3932 sa=sin(ang);
3933 s.x=r0*ca;
3934 s.y=r0*sa;
3935 s.z=z0;
3936 BodyToFrameTAM(&s,&tam,&f);
3937 FrameToBodyTAM(&f,vu,&s);
3938 NormPoint3d(&s,&s,&d);
3939 ToPolar3d(&s,&a);
3940 if(a.x<az->x) az->x=a.x;
3941 if(a.x>az->y) az->y=a.x;
3942 if(a.y<el->x) el->x=a.y;
3943 if(a.y>el->y) el->y=a.y;
3944 s.x=r1*ca;
3945 s.y=r1*sa;
3946 s.z=z1;
3947 BodyToFrameTAM(&s,&tam,&f);
3948 FrameToBodyTAM(&f,vu,&s);
3949 NormPoint3d(&s,&s,&d);
3950 ToPolar3d(&s,&a);
3951 if(a.x<az->x) az->x=a.x;
3952 if(a.x>az->y) az->y=a.x;
3953 if(a.y<el->x) el->x=a.y;
3954 if(a.y>el->y) el->y=a.y;
3955 }
3956}
3957
3958/********************************************************************
3959 CThreeD.cp
3960********************************************************************/
3961
3962void CThreeD::Draw()
3963{
3964 Line3d s;
3965 double aEye,dEye,rEye,xTgt,yTgt,zTgt;
3966 double w,h,d,f,ang1,ang2;
3967 Point2d p;
3968 Point3d x,y,z;
3969 double theta;
3970 double phi;
3971 double psi;
3972 double azi,inc;
3973
3974 // outstr("at the start of CThreeD::Draw()\n");
3975 // flush();
3976 AlignCoord();
3977 VuControl();
3978 // outstr("at the exit of CThreeD::VuControl()\n");
3979 // flush();
3980
3981 aEye=Azi;
3982 dEye=Ele;
3983 rEye=Dst;
3984
3985 xTgt=Tgx;
3986 yTgt=Tgy;
3987 zTgt=Tgz;
3988
3989 if(!contAll){
3990 f=rEye*cos(d2r*dEye);
3991 w=f*cos(d2r*aEye)+xTgt;
3992 h=f*sin(d2r*aEye)+yTgt;
3993 d=rEye*sin(d2r*dEye)+zTgt;
3994 rEye=sqrt(w*w+h*h+d*d);
3995 aEye=atan2(h,w)/d2r;
3996 dEye=atan2(d,sqrt(w*w+h*h))/d2r;
3997 }
3998
3999 s.sx=rEye*cos(d2r*dEye)*cos(d2r*aEye);
4000 s.sy=rEye*cos(d2r*dEye)*sin(d2r*aEye);
4001 s.sz=rEye*sin(d2r*dEye);
4002
4003 s.ex=xTgt;
4004 s.ey=yTgt;
4005 s.ez=zTgt;
4006 MakeThePolar(&s,&thePolar); // Perspective
4007
4008 // convert to double and prepare for conversion to polar coordinate
4009 w=0.5*Wid;
4010 h=0.5*Hei;
4011 d=sqrt(w*w+h*h);
4012 f=Foc;
4013 focus=Foc;
4014
4015 // outstr("at the start of processing the sides\n");
4016 // flush();
4017
4018 // start with right side of the frame
4019 azi=atan2(w,-h); // azimuth
4020 inc=atan2(d,f); // inclination
4021 RtBm.x=azi;
4022 RtBm.y=1.5707963267-inc;
4023 ClipAngleMatrix(azi,inc,&frameRt,&angFrameRt);
4024
4025 // top side of the frame
4026 azi=atan2(w,h); // azimuth
4027 inc=atan2(d,f); // inclination
4028 TpRt.x=azi;
4029 TpRt.y=1.5707963267-inc;
4030 ClipAngleMatrix(azi,inc,&frameTp,&angFrameTp);
4031
4032 // left side of the frame
4033 azi=atan2(-w,h); // azimuth
4034 inc=atan2(d,f); // inclination
4035 LtTp.x=azi;
4036 LtTp.y=1.5707963267-inc;
4037 ClipAngleMatrix(azi,inc,&frameLt,&angFrameLt);
4038
4039 // bottom side of the frame
4040 azi=atan2(-w,-h); // azimuth
4041 inc=atan2(d,f); // inclination
4042 BmLt.x=azi;
4043 BmLt.y=1.5707963267-inc;
4044 ClipAngleMatrix(azi,inc,&frameBm,&angFrameBm);
4045
4046 // outstr("at the start of MemorySetUp()\n");
4047 // flush();
4048
4049 MemorySetUp();
4050
4051 // outstr("at the start of MakePart2()\n");
4052 // flush();
4053
4054 MakePart2(); // part2 is surfaces to be used in the far-near comp
4055
4056 // outstr("at the start of NearFarRelations()\n");
4057 // flush();
4058
4059 NearFarRelations();
4060
4061 // outstr("at the start of DrawView()\n");
4062 // flush();
4063
4064 DrawView();
4065
4066}
4067
4068/******************************************************************************
4069 AlignCoord
4070
4071 calculate general coordinate
4072******************************************************************************/
4073void CThreeD::AlignCoord(void)
4074{
4075 CMPpart *p1;
4076 CMPpart *p2;
4077 TransAngleEuler a,b,c;
4078 int i;
4079
4080 for(i=0;i<npart;i++){
4081 p1=part[i];
4082 a.x=p1->x0;
4083 a.y=p1->y0;
4084 a.z=p1->z0;
4085 a.theta=0.01745329252*p1->theta0;
4086 a.phi=0.01745329252*p1->phi0;
4087 a.psi=0.01745329252*p1->psi0;
4088 if(i==0) c=a;
4089 else{
4090 p2=part[p1->refPart];
4091 b.x=p2->xg;
4092 b.y=p2->yg;
4093 b.z=p2->zg;
4094 b.theta=p2->thetag;
4095 b.phi=p2->phig;
4096 b.psi=p2->psig;
4097 MoveTwiceTAE(&b,&a,&c);
4098 }
4099 p1->xg=c.x;
4100 p1->yg=c.y;
4101 p1->zg=c.z;
4102 p1->thetag=c.theta;
4103 p1->phig=c.phi;
4104 p1->psig=c.psi;
4105 }
4106}
4107
4108
4109void CThreeD::VuControl (void)
4110
4111{
4112 int i;
4113
4114 if(contAll){
4115 TransAngleEuler tae;
4116 TransAngleMatrix tam;
4117 Point2d az,azm,el,elm;
4118 double a,b,e,d;
4119 int i,j;
4120 CMPpart *p;
4121 Point3d los,tgt;
4122
4123 a=0.017453293*Azi;
4124 e=0.017453293*Ele;
4125 d=Dst;
4126 tae.x=d*cos(e)*cos(a);
4127 tae.y=d*cos(e)*sin(a);
4128 tae.z=d*sin(e);
4129 if(e>=0.){
4130 tae.theta=e;
4131 tae.phi=a-1.5707963267;
4132 tae.psi=-1.5707963267;
4133 }
4134 else{
4135 tae.theta=-e;
4136 tae.phi=a+1.5707963267;
4137 tae.psi=1.5707963267;
4138 }
4139 TAEtoTAM(&tae,&tam);
4140 azm.x=10.;
4141 azm.y=-10.;
4142 elm.x=10.;
4143 elm.y=-10.;
4144 // outstr("in CThreeD::VuControl()\n");
4145 // sprintf(msg,"npart=%i\n",npart);
4146 // outstr(msg);
4147 // flush();
4148 for(i=0;i<npart;i++){
4149 p=part[i];
4150 if(p->showDrawing){
4151 switch(p->kind){
4152 case MPlines:
4153 reinterpret_cast<CMPlines *>(p)->FieldOfView(&tam,&az,&el);
4154 break;
4155 case MPplane:
4156 reinterpret_cast<CMPplane *>(p)->FieldOfView(&tam,&az,&el);
4157 break;
4158 case MPbox:
4159 reinterpret_cast<CMPbox *>(p)->FieldOfView(&tam,&az,&el);
4160 break;
4161 case MPsphere:
4162 reinterpret_cast<CMPsphere *>(p)->FieldOfView(&tam,&az,&el);
4163 break;
4164 case MPbodyOfRot:
4165 reinterpret_cast<CMPbodyOfRot *>(p)->FieldOfView(&tam,&az,&el);
4166 break;
4167 default:
4168 az.x=10.;
4169 az.y=-10.;
4170 el.x=10.;
4171 el.y=-10.;
4172 break;
4173 }
4174 if(az.x<azm.x) azm.x=az.x;
4175 if(az.y>azm.y) azm.y=az.y;
4176 if(el.x<elm.x) elm.x=el.x;
4177 if(el.y>elm.y) elm.y=el.y;
4178 }
4179 }
4180 los.x=d*cos(0.5*(elm.x+elm.y))*cos(0.5*(azm.x+azm.y));
4181 los.y=d*cos(0.5*(elm.x+elm.y))*sin(0.5*(azm.x+azm.y));
4182 los.z=d*sin(0.5*(elm.x+elm.y));
4183 BodyToFrameTAM(&los,&tam,&tgt);
4184 Tgx=tgt.x;
4185 Tgy=tgt.y;
4186 Tgz=tgt.z;
4187 a=Wid/2-5;
4188 a=a/tan(0.5*(azm.y-azm.x));
4189 b=Hei/2-5;
4190 b=b/tan(0.5*(elm.y-elm.x));
4191 Foc=(int)((a<b)? a:b);
4192 }
4193}
4194
4195
4196/********************************************************************
4197 MemorySetUp
4198
4199 Analize itsList and allocate memory for following handles
4200
4201 part1part obtain the 'part1' number for a given 'part' number
4202 partpart1 obtain the 'part' number for a given 'part1' number
4203 part2part1 obtain the 'part2' number for a given 'part1' number
4204 part2 array of surfaces indexed by 'part2' number
4205 local array of surface points referred by part2
4206 viewed local converted to polar coordinate
4207 sorted part 1 sorted to the farthest order
4208
4209 part1part and partpart1 are defined in this routine
4210********************************************************************/
4211void CThreeD::MemorySetUp( void)
4212{
4213 int i,j,k,l,m;
4214 int ptesti;
4215 CMPpart *p;
4216 CMPpart *p0; // the first CMPplane of the consecutive group planes
4217 CMPpart *p1;
4218
4219
4220 part1part=new StartPTE[npart]; // (StartPTE *)NewPtr(sizeof(StartPTE)*npart);
4221
4222 npart1=0; // to be incremented by j
4223 npart2m=0; // to be incremented by k
4224 npoinm=0; // to be incremented by l
4225
4226 p0=NULL;
4227
4228 for(i=0;i<npart;i++){
4229 p=part[i]; // itsList->NthItem(i+1)
4230 if(p->showDrawing){
4231 switch(p->kind){
4232 case MPpoint:
4233 j=0; // points are not drawn
4234 k=0; // no near-far comparison
4235 l=0; // no surface points
4236 break;
4237 case MPlines:
4238 j=reinterpret_cast<CMPlines *>(p)->npoint;
4239 j--; // number of segments is number of points minus one
4240 k=j; // one line segment is one generalized surface
4241 l=j*2; // each generalized surface has two points
4242 break;
4243 case MPplane:
4244 j=1; // a plane is a drawing unit
4245 k=1; // one unit in the far-near comparison
4246 l=reinterpret_cast<CMPplane *>(p)->npoint; // points defining the plane
4247
4248 if(reinterpret_cast<CMPplane *>(p)->isGroup){
4249 if(p0==NULL) p0=p;
4250 p1=part[i+1];
4251 if(p1->kind==MPplane){
4252 if(reinterpret_cast<CMPplane *>(p1)->isGroup)
4253 reinterpret_cast<CMPplane *>(p)->otherHalf=reinterpret_cast<CMPplane *>(p1);
4254 else{
4255 reinterpret_cast<CMPplane *>(p)->otherHalf=reinterpret_cast<CMPplane *>(p0);
4256 p0=NULL;
4257 }
4258 }
4259 else{
4260 reinterpret_cast<CMPplane *>(p)->otherHalf=reinterpret_cast<CMPplane *>(p0);
4261 p0=NULL;
4262 }
4263 }
4264 if(reinterpret_cast<CMPplane *>(p)->isPolyhed)
4265 if(!reinterpret_cast<CMPplane *>(p)->FacingViewer((Point3d *)&thePolar)){
4266 j=0;
4267 k=0;
4268 l=0;
4269 }
4270 break;
4271 case MPbox:
4272 j=1; // a box is a drawing unit
4273 k=3; // 3 surfaces max
4274 l=12; // 4 points * 3 surfaces
4275 break;
4276 case MPsphere:
4277 j=1; // a sphere is a drawing unit
4278 k=1; // represented by a surface
4279 l=12; // represented by a 12 sided polygon
4280 break;
4281 case MPbodyOfRot:
4282 if(reinterpret_cast<CMPbodyOfRot *>(p)->isTube){
4283 j=2; // near and far sides
4284 k=12; // 12 surfaces combined, 6 surfaces each
4285 l=48; // 4 points * 12 surfaces
4286 }
4287 else{
4288 j=1; // represented by a drawing unit
4289 k=13; // 12 surfaces + top
4290 l=60; // 4 points * 12 surfaces + 12 points
4291 }
4292 break;
4293 }
4294 }
4295 else{ // if showDrawing is false skip data
4296 j=0;
4297 k=0;
4298 l=0;
4299 }
4300 part1part[i].start=npart1;
4301 npart1+=j;
4302 part1part[i].pte=npart1;
4303 npart2m+=k;
4304 npoinm+=l;
4305 }
4306 npoinm+=5; // alloance for increase by clipping
4307 partpart1=new StartPTE[npart1]; // (StartPTE *)NewPtr(sizeof(StartPTE)*npart1);
4308 for(i=0;i<npart;i++){
4309 j=part1part[i].start;
4310 k=part1part[i].pte;
4311 m=0;
4312 for(l=j;l<k;l++){
4313 partpart1[l].start=i;
4314 partpart1[l].pte=m++;
4315 }
4316 }
4317
4318 part2part1=new StartPTE[npart1]; // (StartPTE *)NewPtr(sizeof(StartPTE)*npart1);
4319 part2=new Surface[npart2m]; // (Surface *)NewPtr(sizeof(Surface)*npart2m);
4320 part1part2=new int[npart2m]; // (int *)NewPtr(sizeof(short)*npart2m);
4321 local=new Point2d[npoinm]; // (Point2d *)NewPtr(sizeof(Point2d)*npoinm);
4322 viewed=new Point2d[npoinm]; //(Point2d *)NewPtr(sizeof(Point2d)*npoinm);
4323 // sorted=(int *)NewPtr(sizeof(short)*npart1);
4324
4325
4326}
4327
4328
4329/********************************************************************
4330 MakePart2
4331********************************************************************/
4332void CThreeD::MakePart2( void)
4333{
4334 int i,j,k,l,m,n;
4335 int ptesti;
4336 CMPpart *p;
4337 double x,y,z,x1,y1,z1;
4338 Surface surf1;
4339 Surface surf2;
4340 Point2d *surfPt1;
4341 Point2d *surfPt2;
4342 Point2d a;
4343 TransAngleMatrix tam;
4344
4345 npart2=0; // number of parts after no 2 decomposition
4346 npoin=0;
4347 surfPt1=new Point2d[50]; // (Point2d *)NewPtr(sizeof(Point2d)*50);
4348 surfPt2=new Point2d[50]; // (Point2d *)NewPtr(sizeof(Point2d)*50);
4349
4350 // outstr("in MakePart2()\n");
4351 // sprintf(msg,"npart=%i\n",npart);
4352 // outstr(msg);
4353 // flush();
4354
4355 for(i=0;i<npart;i++){
4356 j=part1part[i].start;
4357 k=part1part[i].pte;
4358 // sprintf(msg,"i=%i j=%i k=%i\n",i,j,k);
4359 // outstr(msg);
4360 // flush();
4361 if(j<k){
4362 p=part[i]; // (CMPpart *)itsList->NthItem(i+1);
4363 switch(p->kind){
4364 case MPpoint:
4365 part2part1[j].start=-1; // part2 not generated
4366 part2part1[j].pte=-1;
4367 break;
4368 case MPlines:
4369 l=reinterpret_cast<CMPlines *>(p)->npoint;
4370 for(m=0;m<l-1;m++){
4371 reinterpret_cast<CMPlines *>(p)->GetSurface(m,&thePolar,&surf1,surfPt1);
4372 if(surf1.start<0){ // degenerate to a point
4373 surf2.start=-1;
4374 surf2.pte=-1;
4375 }
4376 else ClipSurface(&surf1,surfPt1,&surf2,surfPt2);
4377 if(surf2.start!=surf2.pte){
4378 part2[npart2]=surf2;
4379 part2[npart2].start=npoin;
4380 local[npoin++]=surfPt2[0];
4381 local[npoin++]=surfPt2[1];
4382 part2[npart2].pte=npoin;
4383 part1part2[npart2]=j;
4384 part2part1[j].start=npart2;
4385 npart2++;
4386 part2part1[j].pte=npart2;
4387 }
4388 else{
4389 part2part1[j].start=-1; // part2 not generated
4390 part2part1[j].pte=-1;
4391 }
4392 j++;
4393 }
4394 break;
4395 case MPplane:
4396 reinterpret_cast<CMPplane *>(p)->GetSurface(&thePolar,&surf1,surfPt1);
4397 ClipSurface(&surf1,surfPt1,&surf2,surfPt2);
4398 if(surf2.start!=surf2.pte){
4399 part2[npart2]=surf2;
4400 part2[npart2].start=npoin;
4401 for(n=surf2.start;n<surf2.pte;n++)
4402 local[npoin++]=surfPt2[n];
4403 part2[npart2].pte=npoin;
4404 part1part2[npart2]=j;
4405 part2part1[j].start=npart2;
4406 npart2++;
4407 part2part1[j].pte=npart2;
4408 }
4409 else{
4410 part2part1[j].start=-1; // part2 not generated
4411 part2part1[j].pte=-1;
4412 }
4413 break;
4414 case MPbox:
4415 // outstr("in MakePart2() at MPbox\n");
4416 // sprintf(msg,"npart2=%i\n",npart2);
4417 // outstr(msg);
4418 // flush();
4419
4420 part1part2[npart2]=j;
4421 part2part1[j].start=npart2;
4422 reinterpret_cast<CMPbox *>(p)->ClearVisSurface();
4423 for(m=0;m<3;m++){
4424
4425 if(reinterpret_cast<CMPbox *>(p)->GetSurface(m,&thePolar,&surf1,surfPt1)){
4426 // outstr("at the exit of GetSurface()\n");
4427 // sprintf(msg,"m=%i surf1.x0=%8.3f surf1.y0=%8.3f surf1.z0=%8.3f \n",m,surf1.x0,surf1.y0,surf1.z0);
4428 // outstr(msg);
4429 // sprintf(msg," surf1.theta0=%8.3f surf1.phi0=%8.3f surf1.psi0=%8.3f \n",surf1.theta0,surf1.phi0,surf1.psi0);
4430 // outstr(msg);
4431 // flush();
4432
4433 ClipSurface(&surf1,surfPt1,&surf2,surfPt2);
4434 // outstr("at the exit of ClipSurface()\n");
4435 // sprintf(msg," surf2.start=%i surf2.pte=%i npart2=%i npoin=%i j=%i\n",surf2.start,surf2.pte,npart2,npoin,j);
4436 // outstr(msg);
4437 // flush();
4438 if(surf2.start!=surf2.pte){
4439 part2[npart2]=surf2;
4440 part2[npart2].start=npoin;
4441 for(n=surf2.start;n<surf2.pte;n++)
4442 local[npoin++]=surfPt2[n];
4443 part2[npart2].pte=npoin;
4444 part1part2[npart2]=j;
4445 npart2++;
4446 }
4447 }
4448
4449 }
4450 if(part2part1[j].start!=npart2)
4451 part2part1[j].pte=npart2;
4452 else{
4453 part2part1[j].start=-1; // part2 not generated
4454 part2part1[j].pte=-1;
4455 }
4456 break;
4457 case MPsphere:
4458 reinterpret_cast<CMPsphere *>(p)->GetSurface(&thePolar,&surf1,surfPt1);
4459 ClipSurface(&surf1,surfPt1,&surf2,surfPt2);
4460 if(surf2.start!=surf2.pte){
4461 part2[npart2]=surf2;
4462 part2[npart2].start=npoin;
4463 for(n=surf2.start;n<surf2.pte;n++)
4464 local[npoin++]=surfPt2[n];
4465 part2[npart2].pte=npoin;
4466 part1part2[npart2]=j;
4467 part2part1[j].start=npart2;
4468 npart2++;
4469 part2part1[j].pte=npart2;
4470 }
4471 else{
4472 part2part1[j].start=-1; // part2 not generated
4473 part2part1[j].pte=-1;
4474 }
4475 break;
4476 case MPbodyOfRot:
4477 part1part2[npart2]=j;
4478 part2part1[j].start=npart2;
4479 reinterpret_cast<CMPbodyOfRot *>(p)->ClearVisSurface();
4480 for(m=0;m<13;m++){
4481 if(reinterpret_cast<CMPbodyOfRot *>(p)->GetSurface(m,true, // near side
4482 &thePolar,&surf1,surfPt1)){
4483 ClipSurface(&surf1,surfPt1,&surf2,surfPt2);
4484 if(surf2.start!=surf2.pte){
4485 part2[npart2]=surf2;
4486 part2[npart2].start=npoin;
4487 for(n=surf2.start;n<surf2.pte;n++)
4488 local[npoin++]=surfPt2[n];
4489 part2[npart2].pte=npoin;
4490 part1part2[npart2]=j;
4491 npart2++;
4492 }
4493 }
4494 }
4495 if(part2part1[j].start!=npart2)
4496 part2part1[j].pte=npart2;
4497 else{
4498 part2part1[j].start=-1; // part2 not generated
4499 part2part1[j].pte=-1;
4500 }
4501 if(reinterpret_cast<CMPbodyOfRot *>(p)->isTube){
4502 j++;
4503 part1part2[npart2]=j;
4504 part2part1[j].start=npart2;
4505 for(m=0;m<10;m++){
4506 if(reinterpret_cast<CMPbodyOfRot *>(p)->GetSurface(m,false,
4507 &thePolar,&surf1,surfPt1)){
4508 ClipSurface(&surf1,surfPt1,&surf2,surfPt2);
4509 if(surf2.start!=surf2.pte){
4510 part2[npart2]=surf2;
4511 part2[npart2].start=npoin;
4512 for(n=surf2.start;n<surf2.pte;n++)
4513 local[npoin++]=surfPt2[n];
4514 part1part2[npart2]=j;
4515 part2[npart2++].pte=npoin;
4516 }
4517 }
4518 }
4519 if(part2part1[j].start!=npart2)
4520 part2part1[j].pte=npart2;
4521 else{
4522 part2part1[j].start=-1; // part2 not generated
4523 part2part1[j].pte=-1;
4524 }
4525 }
4526 break;
4527 }
4528 }
4529 }
4530 delete[] surfPt1;
4531 delete[] surfPt2;
4532 // convert local to viewed
4533 for(i=0;i<npart2;i++){
4534 surf1=part2[i];
4535 SurfthePolarTAM((TransAngleEuler *)&surf1,&thePolar,&tam);
4536 for(j=surf1.start;j<surf1.pte;j++){
4537 SurfToPolar(&local[j],&tam,&a);
4538 PolarToGnomo(&a,&viewed[j]);
4539 }
4540 }
4541 delete[] local; // allocated at MemorySetUp()
4542}
4543
4544
4545/********************************************************************
4546 NearFarRelations
4547
4548 establish near-far relationship
4549********************************************************************/
4550
4551void CThreeD::NearFarRelations(void)
4552{
4553 int i,j,k,l,m;
4554 CMPpart *p;
4555 Surface s1,s2;
4556 bool s1isLine;
4557 bool doit;
4558 Point2d p0,p1;
4559 TransAngleMatrix m1,m2;
4560 double d;
4561 char ss[100];
4562 int kind1,kind2;
4563 Point3d center;
4564 Point3d s1z,s2z; // normal to the surface vector in view coordinate
4565 double s1m,s2m; // multiplier in the expression calculating z
4566 long prog,progm;
4567 int progs;
4568 CFarNear *farN;
4569
4570 farN=new CFarNear(npart1);
4571
4572 // enter far-near relationship of tubes
4573 for(i=0;i<npart;i++){
4574 p=part[i];
4575 if(p->showDrawing&&p->kind==MPbodyOfRot){
4576 if(reinterpret_cast<CMPbodyOfRot *>(p)->isTube){
4577 j=part1part[i].start;
4578 farN->AddFarNear(j+1,j);
4579 // outstr("in NearFarRelations()\n");
4580 // sprintf(msg,"i=%i j=%i\n",i,j);
4581 // outstr(msg);
4582 // flush();
4583
4584 }
4585 }
4586 }
4587
4588 // if a line coinsides a side of a plane, put the line in farther position
4589 for(i=0;i<npart;i++){
4590 p=part[i];
4591 if(p->showDrawing&&p->kind==MPplane){
4592 double s1z[3],s20[3];
4593 k=part1part[i].start;
4594 s1=part2[part2part1[k].start];
4595 s1z[0]=sin(s1.theta0)*sin(s1.phi0); // axis normal to the plane
4596 s1z[1]=-sin(s1.theta0)*cos(s1.phi0);
4597 s1z[2]=cos(s1.theta0);
4598 for(j=0;j<npart2;j++){
4599 s2=part2[j];
4600 if(s2.pte-s2.start==2){ // line
4601 s20[0]=s2.x0-s1.x0;
4602 s20[1]=s2.y0-s1.y0;
4603 s20[2]=s2.z0-s1.z0;
4604 d=s1z[0]*s20[0]+s1z[0]*s20[0]+s1z[0]*s20[0];
4605 if(fabs(d)<1.e-10){
4606 m=0;
4607 for(l=s1.start;l<s1.pte;l++){
4608 if(fabs(viewed[i].x-viewed[s2.start].x)<1.e-5
4609 &&fabs(viewed[i].y-viewed[s2.start].y)<1.e-5) m++;
4610 if(fabs(viewed[i].x-viewed[s2.start+1].x)<1.e-5
4611 &&fabs(viewed[i].y-viewed[s2.start+1].y)<1.e-5) m++;
4612 }
4613 if(m==2){
4614 farN->AddFarNear(part1part2[j],k); // line plane
4615 }
4616 }
4617 }
4618 }
4619 }
4620 }
4621
4622 // outstr("in NearFarRelations()\n");
4623 // sprintf(msg,"npart2=%i\n",npart2);
4624 // outstr(msg);
4625 // flush();
4626
4627 // compare each surfaces in the combination triangle
4628 prog=0; // for progress status (long)
4629 progm=npart2;
4630 progm=progm*(progm-1)/2;
4631 for(i=0;i<npart2-1;i++){
4632 k=part1part2[i];
4633 s1=part2[i];
4634 if(s1.pte-s1.start==2) s1isLine=true;
4635 else s1isLine=false;
4636
4637 // sprintf(msg,"i=%i k=%i s1.start=%i s1.pte=%i partpart1[k].start=%i\n",i,k,s1.start,s1.pte,partpart1[k].start);
4638 // outstr(msg);
4639 // flush();
4640
4641 p=part[partpart1[k].start];
4642 kind1=p->kind;
4643 if(kind1==MPbox) reinterpret_cast<CMPbox *>(p)->GetCenter(¢er);
4644 if(kind1==MPbodyOfRot) reinterpret_cast<CMPbodyOfRot *>(p)->GetCenter(¢er);
4645 SurfthePolarTAM((TransAngleEuler *)&s1,&thePolar,&m1);
4646 if(kind1==MPplane){ // used when both parts are planes
4647 Point3d p3;
4648 p3.x=0.;
4649 p3.y=0.;
4650 p3.z=1.; // a point along body z axis
4651 BodyToFrameTAM(&p3,&m1,&s1z);
4652 s1z.x-=m1.x;
4653 s1z.y-=m1.y;
4654 s1z.z-=m1.z;
4655 s1m=s1z.x*m1.x+s1z.y*m1.y+s1z.z*m1.z; // multiplier
4656 }
4657 for(j=i+1;j<npart2;j++){
4658 prog++;
4659 l=part1part2[j];
4660 s2=part2[j];
4661 doit=true;
4662 if(k==l) doit=false;
4663 else if(s1isLine&&(s2.pte-s2.start==2)) doit=false; // dont do it for lines
4664 else doit=!(farN->InferFarNear(k,l,&m));
4665 // sprintf(msg,"j=%i k=%i l=%i doit=%i\n",j,k,l,doit);
4666 // outstr(msg);
4667 //flush();
4668
4669 if(doit){
4670 // sprintf(msg,"part2[i].start=%i part2[i].pte=%i part2[j].start=%i part2[j].pte=%i\n",
4671 // part2[i].start,part2[i].pte,part2[j].start,part2[j].pte);
4672 // outstr(msg);
4673 // flush();
4674
4675 if(Overlap(viewed,(StartPTE *)&part2[i].start,
4676 (StartPTE *)&part2[j].start,&p0)){
4677 // outstr("Overlap\n");
4678 // flush();
4679 SurfthePolarTAM((TransAngleEuler *)&s2,&thePolar,&m2);
4680 GnomoToPolar(&p0,&p1);
4681 FarNearDist(&m1,&m2,&p1,&d);
4682 if(d>0.001){
4683 farN->AddFarNear(k,l);
4684 doit=false;
4685 }
4686 if(d<-0.001){
4687 farN->AddFarNear(l,k);
4688 doit=false;
4689 }
4690 if(doit){
4691 if(kind1==MPplane||kind1==MPbox||kind1==MPbodyOfRot){
4692 p=part[partpart1[l].start];
4693 kind2=p->kind;
4694 if(kind1==MPplane){
4695 if(kind2==MPplane){ // process the planes with common side
4696 Point3d p3;
4697 double z1,z2;
4698 int m,n10,n1p,n1n,n20,n2p,n2n; // zero, plus or negative
4699 p3.x=0.;
4700 p3.y=0.;
4701 p3.z=1.;
4702 BodyToFrameTAM(&p3,&m2,&s2z);
4703 s2z.x-=m2.x;
4704 s2z.y-=m2.y;
4705 s2z.z-=m2.z;
4706 s2m=s2z.x*m2.x+s2z.y*m2.y+s2z.z*m2.z;
4707 n10=0;
4708 n1p=0;
4709 n1n=0;
4710 n20=0;
4711 n2p=0;
4712 n2n=0;
4713 if(fabs(s1z.x*s2z.x+s1z.y*s2z.y+s1z.z*s2z.z)<0.98){ // not parallel
4714 for(m=s1.start;m<s1.pte;m++){
4715 z1=s1m/(s1z.x*viewed[m].y+s1z.y*viewed[m].x+s1z.z);
4716 z2=s2m/(s2z.x*viewed[m].y+s2z.y*viewed[m].x+s2z.z);
4717 if(z2>0.){
4718 if(fabs(z1-z2)<0.001) n10++;
4719 else if(z1>z2) n1p++;
4720 else n1n++;
4721 }
4722 }
4723 for(m=s2.start;m<s2.pte;m++){
4724 z1=s1m/(s1z.x*viewed[m].y+s1z.y*viewed[m].x+s1z.z);
4725 z2=s2m/(s2z.x*viewed[m].y+s2z.y*viewed[m].x+s2z.z);
4726 if(z1>0.){
4727 if(fabs(z1-z2)<0.001) n20++;
4728 else if(z1>z2) n2p++;
4729 else n2n++;
4730 }
4731 }
4732 if(n10>1&&n20>1){
4733 if(n1p>0&&n2p>0){ // s1 is farther
4734 farN->AddFarNear(k,l);
4735 doit=false;
4736 }
4737 if(n1n>0&&n2n>0){ // s2 is farther
4738 farN->AddFarNear(l,k);
4739 doit=false;
4740 }
4741 } // if(n10>1&&n20>1){
4742 } // if(fabs(s1z.x*s2z.x+s1z.y*s2z.y+s1z.z*s2z.z)<0.98){
4743 } // if(kind2==MPplane){
4744 else if(kind2==MPbox||kind2==MPbodyOfRot){
4745 if(kind2==MPbox) reinterpret_cast<CMPbox *>(p)->GetCenter(¢er);
4746 if(kind2==MPbodyOfRot) reinterpret_cast<CMPbodyOfRot *>(p)->GetCenter(¢er);
4747 if(SameSideOfSurface(&s1,¢er,(Point3d *)&thePolar)){
4748 farN->AddFarNear(k,l);
4749 doit=false;
4750 }
4751 else{
4752 farN->AddFarNear(l,k);
4753 doit=false;
4754 }
4755 } // else if(kind2==MPbox||kind2==MPbodyOfRot){
4756 } // if(kind1==MPplane){
4757 else{ // kind1==MPbox or MPbodyOfRot
4758 if(kind2==MPplane){
4759 if(SameSideOfSurface(&s2,¢er,(Point3d *)&thePolar)){
4760 farN->AddFarNear(l,k);
4761 doit=false;
4762 }
4763 else{
4764 farN->AddFarNear(k,l);
4765 doit=false;
4766 }
4767 } // if(kind2==MPplane){
4768 } // else
4769 } // if(kind1==MPplane||kind1==MPbox||kind1==MPbodyOfRot){
4770 } // if(doit){
4771 /* if(doit){
4772
4773 GnomoToPolar(&p0,&p1);
4774 FarNearDist(&m1,&m2,&p1,&d);
4775 if(d>0.){
4776 farN->AddFarNear(k,l);
4777 }
4778 if(d<0.){
4779 farN->AddFarNear(l,k);
4780 }
4781 }*/
4782 } // if(Overlap(viewed,(StartPTE *)&part2[i].start,(StartPTE *)&part2[j].start,&p0)){
4783 } // if(doit){
4784
4785 } // for(j=i+1;j<npart2;j++){
4786 } // for(i=0;i<npart2-1;i++){
4787
4788 delete[] part2part1; // allocated in MemorySetUp
4789 delete[] part1part2; // allocated in MemorySetUp
4790 delete[] part2; // allocated in MemorySetUp
4791 delete[] viewed; // allocated in MemorySetUp
4792 delete[] part1part; // allocated in MemorySetUp at the beginning
4793
4794 sorted=new int[npart1]; // (int *)NewPtr(sizeof(short)*npart1);
4795
4796 m=farN->SetSorted(sorted);
4797 if(m!=0){
4798 // sprintf(ss,"\rloop found in farNear m=%i",m);
4799 // Monitor(90,(Ptr)ss,true);
4800 }
4801}
4802
4803
4804/********************************************************************
4805 DrawView
4806
4807 draw in the view
4808********************************************************************/
4809void CThreeD::DrawView(void)
4810{
4811 int i,j,k,l,m,n,o,q;
4812 StartPTE pt;
4813 CMPpart *p;
4814 Line2d *seg;
4815 StartPTE blk;
4816 StartPTE cont;
4817 double x,y,r,ang1,ang2;
4818 int ix,iy;
4819 double maxang;
4820 int red,green,blue;
4821 bool finished[200];
4822 bool isLines;
4823 double a[6];
4824 int b[4];
4825 FILE *fp;
4826 char *buff;
4827 int repeat;
4828 FILE *fp1;
4829 char c;
4830 Point2d p2d1,p2d2;
4831 Point3d p3d1,p3d2,p3d3;
4832 AngleMatrix am1;
4833 double w,h,d,f,azi,inc;
4834
4835 maxang=atan(10000./focus);
4836
4837 seg=new Line2d[200]; // (Line2d *)NewPtr(sizeof(Line2d)*200);
4838
4839 if(postscript){
4840 fp1=fopen("/tmp/gcs3DPOST","w");
4841 if(postscript==1){
4842 fputs("%!PS-Adobe-3.0 EPSF-3.0\n%%BoundingBox: 0 0 ",fp1);
4843 sprintf(msg,"%i %i\n",2*Wid,2*Hei);
4844 fputs(msg,fp1);
4845 }
4846 else if(postscript==2){
4847 fputs("<html>\n<head>\n<script type=\"application/x-javascript\">\nfunction draw() {\n"
4848 "var canvas = document.getElementById(\"canvas\");\nvar ctx = canvas.getContext(\"2d\");\n",fp1);
4849 }
4850 else{
4851 fputs("<?xml version=\"1.0\"?>\n<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n"
4852 "\"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n<svg xmlns=\"http://www.w3.org/2000/svg\"\n"
4853 "xmlns:xlink=\"http://www.w3.org/1999/xlink\"\n",fp1);
4854 sprintf(msg," width=\"%i\" height=\"%i\">\n",2*Wid,2*Hei);
4855 fputs(msg,fp1);
4856 fputs("<script>\nvar color;\nfunction On(a){\na.target.setAttribute('fill','#ff0000');\n"
4857 "}\nfunction Off(a){\na.target.setAttribute('fill',color);\n}\n</script>\n",fp1);
4858 }
4859 }
4860 else if(omni==0){
4861 MakeView(2*Wid,2*Hei);
4862 SetColor(100,100,100);
4863 OpenPoly();
4864 MoveTo(0,0);
4865 LineTo(0,2*Hei);
4866 LineTo(2*Wid,2*Hei);
4867 LineTo(2*Wid,0);
4868 ClosePoly();
4869 }
4870 else{
4871 MakeView(602,301);
4872 SetColor(100,100,100);
4873 OpenPoly();
4874 MoveTo(0,0);
4875 LineTo(0,301);
4876 LineTo(602,301);
4877 LineTo(602,0);
4878 ClosePoly();
4879 SetColor(100,60,30);
4880 MoveTo(150+150,150);
4881 for(x=0.;x<6.3;x+=0.03141592)
4882 LineTo((int)(150+150.*cos(x)),(int)(150-150.*sin(x)));
4883 MoveTo(451+150,150);
4884 for(x=0.;x<6.3;x+=0.03141592)
4885 LineTo((int)(451+150.*cos(x)),(int)(150-150.*sin(x)));
4886 SetColor(100,0,0);
4887 r=150.*(1.570796327-RtBm.y)/1.570796327;
4888 MoveTo((int)(150+r*sin(RtBm.x)),(int)(150-r*cos(RtBm.x)));
4889 for(ang1=0.;ang1<angFrameRt;ang1+=0.03){
4890 p2d2.x=ang1;
4891 p2d2.y=0.;
4892 BodyToFrameAMP(&p2d2,&frameRt,&p2d1);
4893 r=150.*(1.570796327-p2d1.y)/1.570796327;
4894 LineTo((int)(150+r*sin(p2d1.x)),(int)(150-r*cos(p2d1.x)));
4895 }
4896 for(ang1=0.;ang1<angFrameTp;ang1+=0.03){
4897 p2d2.x=ang1;
4898 p2d2.y=0.;
4899 BodyToFrameAMP(&p2d2,&frameTp,&p2d1);
4900 r=150.*(1.570796327-p2d1.y)/1.570796327;
4901 LineTo((int)(150+r*sin(p2d1.x)),(int)(150-r*cos(p2d1.x)));
4902 }
4903 for(ang1=0.;ang1<angFrameLt;ang1+=0.03){
4904 p2d2.x=ang1;
4905 p2d2.y=0.;
4906 BodyToFrameAMP(&p2d2,&frameLt,&p2d1);
4907 r=150.*(1.570796327-p2d1.y)/1.570796327;
4908 LineTo((int)(150+r*sin(p2d1.x)),(int)(150-r*cos(p2d1.x)));
4909 }
4910 for(ang1=0.;ang1<angFrameBm;ang1+=0.03){
4911 p2d2.x=ang1;
4912 p2d2.y=0.;
4913 BodyToFrameAMP(&p2d2,&frameBm,&p2d1);
4914 r=150.*(1.570796327-p2d1.y)/1.570796327;
4915 LineTo((int)(150+r*sin(p2d1.x)),(int)(150-r*cos(p2d1.x)));
4916 }
4917 r=150.*(1.570796327-RtBm.y)/1.570796327;
4918 LineTo((int)(150+r*sin(RtBm.x)),(int)(150-r*cos(RtBm.x)));
4919 SetColor(0,0,0);
4920 }
4921
4922 if((omni&1)==0){ /* clipping - parenteses are required */
4923 // convert to double and prepare for conversion to polar coordinate
4924 w=0.55*Wid; /* make clipping rect slightly larger than the real frame */
4925 h=0.55*Hei;
4926 d=sqrt(w*w+h*h);
4927 f=Foc;
4928
4929 // start with right side of the frame
4930 azi=atan2(w,-h); // azimuth
4931 inc=atan2(d,f); // inclination
4932 RtBm.x=azi;
4933 RtBm.y=1.5707963267-inc;
4934 ClipAngleMatrix(azi,inc,&frameRt,&angFrameRt);
4935
4936 // top side of the frame
4937 azi=atan2(w,h); // azimuth
4938 inc=atan2(d,f); // inclination
4939 TpRt.x=azi;
4940 TpRt.y=1.5707963267-inc;
4941 ClipAngleMatrix(azi,inc,&frameTp,&angFrameTp);
4942
4943 // left side of the frame
4944 azi=atan2(-w,h); // azimuth
4945 inc=atan2(d,f); // inclination
4946 LtTp.x=azi;
4947 LtTp.y=1.5707963267-inc;
4948 ClipAngleMatrix(azi,inc,&frameLt,&angFrameLt);
4949
4950 // bottom side of the frame
4951 azi=atan2(-w,-h); // azimuth
4952 inc=atan2(d,f); // inclination
4953 BmLt.x=azi;
4954 BmLt.y=1.5707963267-inc;
4955 ClipAngleMatrix(azi,inc,&frameBm,&angFrameBm);
4956 }
4957
4958 for(i=0;i<npart1;i++){
4959 j=sorted[i]; // part 1 index
4960 pt=partpart1[j]; // part index
4961 k=pt.start;
4962 p=part[k]; // (CMPpart *)itsList->NthItem(k+1);
4963 isLines=false;
4964 repeat=1;
4965 switch (p->kind){
4966 case MPlines:
4967 reinterpret_cast<CMPlines *>(p)->GetDrawData(pt.pte,&thePolar,seg,&blk,&cont);
4968 isLines=true;
4969 break;
4970 case MPplane:
4971 reinterpret_cast<CMPplane *>(p)->GetDrawData(&thePolar,seg,&blk,&cont);
4972 break;
4973 case MPbox:
4974 l=reinterpret_cast<CMPbox *>(p)->objid;
4975 if(l==0) reinterpret_cast<CMPbox *>(p)->GetDrawData(&thePolar,seg,&blk,&cont);
4976 else{
4977 a[0]=Azi;
4978 a[1]=Ele;
4979 a[2]=Dst;
4980 a[3]=Tgx;
4981 a[4]=Tgy;
4982 a[5]=Tgz;
4983 b[0]=Wid;
4984 b[1]=Hei;
4985 b[2]=Foc;
4986 b[3]=contAll;
4987 setting(a,b);
4988 a[0]=p->xg;
4989 a[1]=p->yg;
4990 a[2]=p->zg;
4991 a[3]=p->thetag*57.29577951;
4992 a[4]=p->phig*57.29577951;
4993 a[5]=p->psig*57.29577951;
4994 drawobj(l,a);
4995 fp=fopen("/tmp/gcs3DDATA","r");
4996 fread(&l,sizeof(int),1,fp); // npart
4997 fread(&l,sizeof(int),1,fp); // npart1
4998 repeat=2;
4999 }
5000 break;
5001 case MPsphere:
5002 reinterpret_cast<CMPsphere *>(p)->GetDrawData(&thePolar,seg,&blk,&cont);
5003 break;
5004 case MPbodyOfRot:
5005 reinterpret_cast<CMPbodyOfRot *>(p)->GetDrawData(pt.pte,&thePolar,seg,&blk,&cont);
5006 break;
5007 }
5008 while(1){
5009 if(repeat==2){
5010 fread(&k,sizeof(int),1,fp); // size of buff
5011 if(k==0) break;
5012 else{
5013 buff=(char *)malloc(sizeof(char)*k);
5014 fread(buff,sizeof(char),k,fp);
5015 p->color=*(int *)buff;
5016 isLines=*(int *)(buff+4);
5017 blk.start=*(int *)(buff+8);
5018 l=blk.pte=*(int *)(buff+12);
5019 cont.start=*(int *)(buff+16);
5020 cont.pte=*(int *)(buff+20);
5021 for(j=0;j<l;j++){
5022 seg[j].sx=*(double *)(buff+24+j*32);
5023 seg[j].sy=*(double *)(buff+24+j*32+8);
5024 seg[j].ex=*(double *)(buff+24+j*32+16);
5025 seg[j].ey=*(double *)(buff+24+j*32+24);
5026 }
5027 }
5028 free(buff);
5029 } /* if(repeat==2){ */
5030 if((omni&1)==0){ /* clipping */
5031 ClipSegments(&blk,&cont,seg);
5032 }
5033 if(omni){
5034 k=cont.start;
5035 l=blk.pte;
5036 for(j=k;j<l;j++){
5037 p3d1.z=sin(seg[j].sy);
5038 p3d1.x=p3d1.z*cos(seg[j].sx);
5039 p3d1.y=p3d1.z*sin(seg[j].sx);
5040 p3d1.z=cos(seg[j].sy);
5041 p3d2.z=sin(seg[j].ey);
5042 p3d2.x=p3d2.z*cos(seg[j].ex);
5043 p3d2.y=p3d2.z*sin(seg[j].ex);
5044 p3d2.z=cos(seg[j].ey);
5045 ang1=acos(p3d1.x*p3d2.x+p3d1.y*p3d2.y+p3d1.z*p3d2.z);
5046 if(ang1!=0.){
5047 VectorProduct(&p3d1,&p3d2,&p3d3);
5048 PointsToAM(&p3d3,&p3d1,zxAxis,&am1);
5049 p2d1.x=seg[j].sx;
5050 p2d1.y=seg[j].sy;
5051 if(p2d1.y<1.570796327){
5052 r=p2d1.y*150/1.570796327;
5053 MoveTo((int)(150+r*sin(p2d1.x)),(int)(150-r*cos(p2d1.x)));
5054 }
5055 else{
5056 r=(3.141592654-p2d1.y)*150/1.570796327;
5057 MoveTo((int)(450-r*sin(p2d1.x)),(int)(150-r*cos(p2d1.x)));
5058 }
5059 for(ang2=0.03;ang2<ang1;ang2+=0.03){
5060 p2d2.x=ang2;
5061 p2d2.y=0.;
5062 BodyToFrameAMP(&p2d2,&am1,&p2d1);
5063 p2d1.y=1.570796327-p2d1.y;
5064 if(p2d1.y<1.570796327){
5065 r=p2d1.y*150/1.570796327;
5066 LineTo((int)(150+r*sin(p2d1.x)),(int)(150-r*cos(p2d1.x)));
5067 }
5068 else{
5069 r=(3.141592654-p2d1.y)*150/1.570796327;
5070 LineTo((int)(450-r*sin(p2d1.x)),(int)(150-r*cos(p2d1.x)));
5071 }
5072 } /* for */
5073 p2d1.x=seg[j].ex;
5074 p2d1.y=seg[j].ey;
5075 if(p2d1.y<1.570796327){
5076 r=p2d1.y*150/1.570796327;
5077 LineTo((int)(150+r*sin(p2d1.x)),(int)(150-r*cos(p2d1.x)));
5078 }
5079 else{
5080 r=(3.141592654-p2d1.y)*150/1.570796327;
5081 LineTo((int)(450-r*sin(p2d1.x)),(int)(150-r*cos(p2d1.x)));
5082 }
5083 } /* if */
5084 } /* for */
5085 if(repeat==2) continue;
5086 break;
5087 } /* if */
5088 k=cont.start; // supposed to cover all segments
5089 l=blk.pte;
5090 for(j=k;j<l;j++){
5091 x=seg[j].sx;
5092 y=seg[j].sy;
5093 if(y>maxang) y=maxang;
5094 y=2*focus*tan(y);
5095 seg[j].sx=y*sin(x);
5096 seg[j].sy=y*cos(x);
5097 x=seg[j].ex;
5098 y=seg[j].ey;
5099 if(y>maxang) y=maxang;
5100 y=2*focus*tan(y);
5101 seg[j].ex=y*sin(x);
5102 seg[j].ey=y*cos(x);
5103 }
5104
5105 if(l-k>1){
5106 switch(p->color){
5107 case 0:
5108 red=100;
5109 green=100;
5110 blue=100;
5111 break;
5112 case 1:
5113 red=75;
5114 green=75;
5115 blue=75;
5116 break;
5117 case 2:
5118 red=100;
5119 green=75;
5120 blue=75;
5121 break;
5122 case 3:
5123 red=100;
5124 green=100;
5125 blue=75;
5126 break;
5127 case 4:
5128 red=75;
5129 green=100;
5130 blue=75;
5131 break;
5132 case 5:
5133 red=75;
5134 green=100;
5135 blue=100;
5136 break;
5137 case 6:
5138 red=75;
5139 green=75;
5140 blue=100;
5141 break;
5142 case 7:
5143 red=100;
5144 green=75;
5145 blue=100;
5146 break;
5147 }
5148 if(postscript){
5149 if(postscript==1) sprintf(msg,"%f %f %f setrgbcolor\n",0.01*red,0.01*green,0.01*blue);
5150 else if(postscript==2) sprintf(msg,"ctx.fillStyle=\"rgb(%i,%i,%i)\";\n",255*red/100,255*green/100,255*blue/100);
5151 else sprintf(msg,"<g style=\"fill: #%02x%02x%02x\">\n<path onmouseover=\"On(evt)\" onmouseout=\"Off(evt)\" d=\"",
5152 255*red/100,255*green/100,255*blue/100);
5153 fputs(msg,fp1);
5154 }
5155 else SetColor(red,green,blue);
5156 // PenSize(1,1);
5157 k=cont.start; // start of contour
5158 l=cont.pte-cont.start; // number of contour segments
5159 q=k;
5160 if(postscript){
5161 if(postscript==1) sprintf(msg,"newpath\n%.2f %.2f moveto\n",seg[q].sx+Wid,Hei+seg[q].sy);
5162 else if(postscript==2) sprintf(msg,"ctx.beginPath();\nctx.moveTo(%.2f,%.2f);\n",seg[q].sx+Wid,Hei-seg[q].sy);
5163 else sprintf(msg,"M%.1f %.1f",seg[q].sx+Wid,Hei-seg[q].sy);
5164 fputs(msg,fp1);
5165 }
5166 else{
5167 OpenPoly();
5168 ix=(int)(seg[q].sx+Wid);
5169 iy=(int)(Hei-seg[q].sy);
5170 MoveTo(ix,iy);
5171 }
5172 for(n=0;n<l;n++) finished[n]=false;
5173 for(m=0;m<l;m++){
5174 if(postscript){
5175 if(postscript==1) sprintf(msg,"%.2f %.2f lineto\n",seg[q].ex+Wid,Hei+seg[q].ey);
5176 else if(postscript==2) sprintf(msg,"ctx.lineTo(%.2f,%.2f);\n",seg[q].ex+Wid,Hei-seg[q].ey);
5177 else sprintf(msg,"L%.1f %.1f",seg[q].ex+Wid,Hei-seg[q].ey);
5178 fputs(msg,fp1);
5179 }
5180 else{
5181 ix=(int)(seg[q].ex+Wid);
5182 iy=(int)(Hei-seg[q].ey);
5183 LineTo(ix,iy);
5184 }
5185 for(n=0;n<l;n++){
5186 if(!finished[n]){
5187 o=k+n;
5188 if(seg[q].ex==seg[o].sx && seg[q].ey==seg[o].sy){
5189 q=o;
5190 finished[n]=true;
5191 break;
5192 }
5193 }
5194 }
5195 }
5196 if(postscript){
5197 if(postscript==1) fputs("closepath\nfill\n",fp1);
5198 else if(postscript==2) fputs("ctx.closePath();\nctx.fill();\n",fp1);
5199 else fputs(" z\"/></g>\n",fp1);
5200 }
5201 else ClosePoly();
5202 } /* if(l-k>1){ */
5203 // PenSize(ContPen,ContPen);
5204 if(postscript){
5205 if(postscript==1) fputs("0 0 0 setrgbcolor\n",fp1);
5206 else if(postscript==2) fputs("ctx.strokeStyle=\"rgb(0,0,0)\";\n",fp1);
5207 else ; // fputs("<g style=\"stroke: #000000; stroke-width:1\"><path d=\"",fp1);
5208 }
5209 else SetColor(0,0,0);
5210 if(isLines){
5211 red=0;
5212 green=0;
5213 blue=0;
5214 switch(p->color){
5215 case 1:
5216 // PenSize(2,2);
5217 break;
5218 case 2:
5219 // PenSize(4,4);
5220 break;
5221 case 3:
5222 // PenSize(8,8);
5223 break;
5224 case 4:
5225 red=100;
5226 green=0;
5227 blue=0;
5228 // PenSize(4,4);
5229 break;
5230 case 5:
5231 red=100;
5232 green=0;
5233 blue=0;
5234 // PenSize(8,8);
5235 break;
5236 case 6:
5237 red=0;
5238 green=0;
5239 blue=100;
5240 // PenSize(4,4);
5241 break;
5242 case 7:
5243 red=0;
5244 green=0;
5245 blue=100;
5246 // PenSize(8,8);
5247 break;
5248 }
5249 if(postscript){
5250 if(postscript==1) sprintf(msg,"%f %f %f setrgbcolor\n",0.01*red,0.01*green,0.01*blue);
5251 else if(postscript==2) sprintf(msg,"ctx.strokeStyle=\"rgb(%i,%i,%i)\";\n",255*red/100,255*green/100,255*blue/100);
5252 else sprintf(msg,"<g style=\"stroke: #%02x%02x%02x; stroke-width:1\"><path d=\"",255*red/100,255*green/100,255*blue/100);
5253 fputs(msg,fp1);
5254 }
5255 else SetColor(red,green,blue);
5256 }
5257 k=blk.start;
5258 l=blk.pte;
5259 for(m=k;m<l;m++){
5260 if(postscript){
5261 if(postscript==1) sprintf(msg,"newpath\n%.2f %.2f moveto\n",seg[m].sx+Wid,Hei+seg[m].sy);
5262 else if(postscript==2) sprintf(msg,"ctx.beginPath();\nctx.moveTo(%.2f,%.2f);\n",seg[m].sx+Wid,Hei-seg[m].sy);
5263 else sprintf(msg,"<g style=\"stroke: #000000\"><path d=\"M%.1f %.1f",seg[m].sx+Wid,Hei-seg[m].sy);
5264 fputs(msg,fp1);
5265 }
5266 else{
5267 ix=(int)(seg[m].sx+Wid);
5268 iy=(int)(Hei-seg[m].sy);
5269 MoveTo(ix,iy);
5270 }
5271 if(postscript){
5272 if(postscript==1) sprintf(msg,"%.2f %.2f lineto\nstroke\n",seg[m].ex+Wid,Hei+seg[m].ey);
5273 else if(postscript==2) sprintf(msg,"ctx.lineTo(%.2f,%.2f);\nctx.stroke();\n",seg[m].ex+Wid,Hei-seg[m].ey);
5274 else sprintf(msg,"L%.1f %.1f\"/></g>\n",seg[m].ex+Wid,Hei-seg[m].ey);
5275 fputs(msg,fp1);
5276 }
5277 else{
5278 ix=(int)(seg[m].ex+Wid);
5279 iy=(int)(Hei-seg[m].ey);
5280 LineTo(ix,iy);
5281 }
5282 }
5283 if(isLines){
5284 // PenSize(1,1);
5285 red=0;
5286 green=0;
5287 blue=0;
5288 if(postscript){
5289 if(postscript==1) sprintf(msg,"%f %f %f setrgbcolor\n",0.01*red,0.01*green,0.01*blue);
5290 else if(postscript==2) sprintf(msg,"ctx.fillstyle=\"rgb(%i,%i,%i)\";",255*red/100,255*green/100,255*blue/100);
5291 /* svg to be implemented */
5292 fputs(msg,fp1);
5293 }
5294 else SetColor(red,green,blue);
5295 }
5296 if(repeat==1) break;
5297 } /* while */
5298 if(repeat==2) fclose(fp);
5299 } /* for(i=0;i<npart1;i++){ */
5300
5301 if(message==1&&postscript==0){
5302 SetColor(0,0,80);
5303 FontInfo(11,0,0);
5304 MoveTo(1,11);
5305 sprintf(msg,"Azi=%8.3f",Azi);
5306 DrawString(msg);
5307 MoveTo(1,11*2);
5308 sprintf(msg,"Ele=%8.3f",Ele);
5309 DrawString(msg);
5310 MoveTo(1,11*3);
5311 sprintf(msg,"Dst=%8.3f",Dst);
5312 DrawString(msg);
5313 MoveTo(1,11*4);
5314 sprintf(msg,"Tgx=%8.3f",Tgx);
5315 DrawString(msg);
5316 MoveTo(1,11*5);
5317 sprintf(msg,"Tgy=%8.3f",Tgy);
5318 DrawString(msg);
5319 MoveTo(1,11*6);
5320 sprintf(msg,"Tgz=%8.3f",Tgz);
5321 DrawString(msg);
5322 MoveTo(1,11*7);
5323 sprintf(msg,"Wid=%5i",Wid*2);
5324 DrawString(msg);
5325 MoveTo(1,11*8);
5326 sprintf(msg,"Hei=%5i",Hei*2);
5327 DrawString(msg);
5328 MoveTo(1,11*9);
5329 sprintf(msg,"Foc=%5i",Foc*2);
5330 DrawString(msg);
5331 }
5332 if(postscript){
5333 if(postscript==1) fputs("showpage\n",fp1);
5334 else if(postscript==2) {
5335 fputs("}\n</script>\n</head>\n<body onload=\"draw()\">\n<canvas id=\"canvas\" ",fp1);
5336 sprintf(msg,"width=\"%i\" height=\"%i\"></canvas>\n</body>\n</html>\n",2*Wid,2*Hei);
5337 fputs(msg,fp1);
5338 }
5339 else fputs("</svg>\n",fp1);
5340 fclose(fp1);
5341 fp1=fopen("/tmp/gcs3DPOST","r");
5342 if(postscript==3){
5343 printf("Content-type: image/svg+xml%c%c",10,10);
5344 while((i=fgetc(fp1))!=EOF){
5345 c=i;
5346 putchar(c);
5347 }
5348 }
5349 else{
5350 while((i=fgetc(fp1))!=EOF){
5351 c=i;
5352 outbyte(1,&c);
5353 }
5354 flush();
5355 }
5356 fclose(fp1);
5357 }
5358 delete[] seg; // DisposPtr((Ptr)seg);
5359 delete[] partpart1; // DisposPtr((Ptr)partpart1);
5360 delete[] sorted; // DisposPtr((Ptr)sorted);
5361
5362
5363}
5364
5365
5366/********************************************************************
5367 ClipSurface
5368
5369 clip a surface by the view frame
5370********************************************************************/
5371void CThreeD::ClipSurface(Surface *surf1,Point2d *surfPt1,
5372 Surface *surf2,Point2d *surfPt2)
5373{
5374 TransAngleMatrix s; // body:surface, frame:thePolar
5375 Point2d p0,p1,p2,p3,p4;
5376 int nx; // number of crossing points
5377 Point2d p[8]; // crossing point
5378 int ix[8]; // start point index
5379 int jx[8]; // status 1 : coming-in, 2 : going-out
5380 int kx[8]; // clip-side 1:Rt, 2:Tp, 3:Lt, 4:Bm
5381 double ang[8]; // angle from the corner
5382 int status; // returned status
5383 double an; // returned angle
5384 bool inside;
5385 int i,j,k,l;
5386
5387 SurfthePolarTAM((TransAngleEuler *)surf1,&thePolar,&s); // Perspective
5388 SurfToPolar(&surfPt1[0],&s,&p0); // (x,y) to (alpha,delta)
5389 p2=p0;
5390 nx=0;
5391 // outstr("in ClipSurface()\n");
5392 // sprintf(msg," surf1->start=%i surf1->pte=%i \n",surf1->start,surf1->pte);
5393 // outstr(msg);
5394 // flush();
5395
5396 for(i=surf1->start;i<surf1->pte;i++){
5397 p1=p2;
5398 if(i+1==surf1->pte){
5399 if(surf1->pte-surf1->start==2) break;
5400 else p2=p0;
5401 }
5402 else SurfToPolar(&surfPt1[i+1],&s,&p2);
5403
5404 status=CrossClipSide(&frameRt,angFrameRt,&p1,&p2,&p3,&an); // Perspective
5405
5406 if(status<3){ // crossing or inside
5407 if(status>0){ // crossing
5408 p[nx]=p3; // crossing point in polar coordinate
5409 ix[nx]=i; // start point index
5410 jx[nx]=status; // status 1 or 2
5411 kx[nx]=1; // clip rect side - right
5412 ang[nx]=an; // distance from right-bottom corner
5413 nx++;
5414 }
5415 status=CrossClipSide(&frameTp,angFrameTp,&p1,&p2,&p3,&an);
5416
5417 }
5418 if(status<3){
5419 if(status>0){
5420 p[nx]=p3;
5421 ix[nx]=i;
5422 jx[nx]=status;
5423 kx[nx]=2;
5424 ang[nx]=an;
5425 nx++;
5426 }
5427 status=CrossClipSide(&frameLt,angFrameLt,&p1,&p2,&p3,&an);
5428
5429 }
5430 if(status<3){
5431 if(status>0){
5432 p[nx]=p3;
5433 ix[nx]=i;
5434 jx[nx]=status;
5435 kx[nx]=3;
5436 ang[nx]=an;
5437 nx++;
5438 }
5439 status=CrossClipSide(&frameBm,angFrameBm,&p1,&p2,&p3,&an);
5440
5441 }
5442 if(status<3){
5443 if(status>0){
5444 p[nx]=p3;
5445 ix[nx]=i;
5446 jx[nx]=status;
5447 kx[nx]=4;
5448 ang[nx]=an;
5449 nx++;
5450 }
5451 }
5452
5453
5454 }
5455 // outstr("in ClipSurface()\n");
5456 // sprintf(msg," nx=%i status=%i surf1->start=%i surf1->pte=%i \n",nx,status,surf1->start,surf1->pte);
5457 // outstr(msg);
5458 // flush();
5459
5460
5461 if(nx==0){
5462 if(status==0){ // the whole surface is inside the clip rect
5463 *surf2=*surf1;
5464 for(i=surf1->start;i<surf1->pte;i++) surfPt2[i]=surfPt1[i];
5465 return;
5466 }
5467 else{
5468 *surf2=*surf1;
5469 surf2->pte=surf2->start; // set past-the-end to start
5470 return;
5471 }
5472 }
5473 // if two crossings occur for a segment, 'enter' comes first and then comes 'leave'
5474 for(i=1;i<nx;i++){
5475 if(ix[i]==ix[i-1]){
5476 if(jx[i]==1){ // enter
5477 p3=p[i];
5478 p[i]=p[i-1];
5479 p[i-1]=p3;
5480 j=jx[i];
5481 jx[i]=jx[i-1];
5482 jx[i-1]=j;
5483 j=kx[i];
5484 kx[i]=kx[i-1];
5485 kx[i-1]=j;
5486 an=ang[i];
5487 ang[i]=ang[i-1];
5488 ang[i-1]=an;
5489 }
5490 }
5491 }
5492
5493 // add corner points after 'leave' and before 'enter'
5494 for(i=0;i<nx;i++){
5495 if(jx[i]==2){ // leave
5496 j=i+1;
5497 if(j==nx) j=0; // next enter point
5498 k=kx[i]; // side of leave point
5499 while(k!=kx[j]){
5500 k++;
5501 if(k==5) k=1;
5502 for(l=nx;l>i+1;l--){ // make room for a corner point
5503 p[l]=p[l-1];
5504 ix[l]=ix[l-1];
5505 jx[l]=jx[l-1];
5506 kx[l]=kx[l-1];
5507 ang[l]=ang[l-1];
5508 }
5509 switch(k){
5510 case 1:
5511 p[i+1]=RtBm;
5512 break;
5513 case 2:
5514 p[i+1]=TpRt;
5515 break;
5516 case 3:
5517 p[i+1]=LtTp;
5518 break;
5519 case 4:
5520 p[i+1]=BmLt;
5521 break;
5522 }
5523 ix[i+1]=ix[i];
5524 jx[i+1]=4; // status is corner point
5525 kx[i+1]=k;
5526 ang[i+1]=0.;
5527 nx++;
5528 i++;
5529 }
5530 }
5531 }
5532
5533 // convert polar coordinate to surface coordinate
5534 for(i=0;i<nx;i++) PolarToSurf(&p[i],&s,&p[i]);
5535
5536 // generate new surface
5537 j=0; // index to crossing points
5538 k=0; // index to output points
5539 for(i=surf1->start;i<surf1->pte;i++){
5540 inside=false;
5541 for(j=0;j<nx;j++){
5542 if(i<=ix[j]&&jx[j]==2) {
5543 inside=true;
5544 break;
5545 }
5546 if(i<=ix[j]&&jx[j]==1) {
5547 inside=false;
5548 break;
5549 }
5550 }
5551 if(j==nx&&jx[0]==2) inside=true;
5552 if(inside){
5553 surfPt2[k]=surfPt1[i];
5554 k++;
5555 }
5556 for(j=0;j<nx;j++){
5557 if(i==ix[j]){
5558 surfPt2[k]=p[j];
5559 k++;
5560 }
5561 }
5562 }
5563 *surf2=*surf1;
5564 surf2->start=0;
5565 surf2->pte=k;
5566
5567}
5568
5569/********************************************************************
5570 ClipSegments
5571
5572 clip line segments by the view frame
5573********************************************************************/
5574void CThreeD::ClipSegments(StartPTE *blk,StartPTE *cont,Line2d *seg)
5575{
5576 Point2d con[200];
5577 Point2d con1[200];
5578 Line2d seg1[200];
5579 Point2d p0,p1,p2,p3,p4;
5580 int nx; // number of crossing points
5581 Point2d p[8]; // crossing point
5582 int ix[8]; // start point index
5583 int jx[8]; // status 1 : coming-in, 2 : going-out
5584 int kx[8]; // clip-side 1:Rt, 2:Tp, 3:Lt, 4:Bm
5585 double ang[8]; // angle from the corner
5586 int status; // returned status
5587 double an; // returned angle
5588 bool inside;
5589 int i,j,k,l,m,n,o,q;
5590 bool finished[200];
5591 double x,y;
5592
5593 /* quickly test if the clipping is required */
5594 k=cont->start;
5595 l=blk->pte;
5596 y=1.5707963267-RtBm.y; /* radius */
5597 for(j=k;j<l;j++){
5598 if(seg[j].sy>y) break;
5599 if(seg[j].ey>y) break;
5600 }
5601 if(j==l) return; /* within close to the boundary */
5602 k=cont->start; // start of contour
5603 l=cont->pte-cont->start; // number of contour segments
5604 if(l!=0){
5605 j=0; /* make the contour in the order */
5606 q=k;
5607 con[j].x=seg[q].sx;
5608 con[j].y=1.570796327-seg[q].sy; /* convert to declination */
5609 j++;
5610 for(n=0;n<l;n++) finished[n]=false;
5611 for(m=0;m<l;m++){
5612 con[j].x=seg[q].ex;
5613 con[j].y=1.570796327-seg[q].ey;
5614 j++;
5615 for(n=0;n<l;n++){
5616 if(!finished[n]){
5617 o=k+n;
5618 if(seg[q].ex==seg[o].sx && seg[q].ey==seg[o].sy){
5619 q=o;
5620 finished[n]=true;
5621 break;
5622 }
5623 }
5624 }
5625 }
5626 if(con[j-1].x==con[0].x && con[j-1].y==con[0].y) j--;
5627 nx=0; /* number of crossing points */
5628 for(i=0;i<j;i++){ /* do for each segment for the contour */
5629 p1=con[i];
5630 if(i==j-1) p2=con[0];
5631 else p2=con[i+1];
5632 status=CrossClipSide(&frameRt,angFrameRt,&p1,&p2,&p3,&an);
5633 if(status<3){ // crossing or inside
5634 if(status>0){ // crossing
5635 p[nx]=p3; // crossing point in polar coordinate
5636 ix[nx]=i; // start point index
5637 jx[nx]=status; // status 1 or 2
5638 kx[nx]=1; // clip rect side - right
5639 ang[nx]=an; // distance from right-bottom corner
5640 nx++;
5641 }
5642 status=CrossClipSide(&frameTp,angFrameTp,&p1,&p2,&p3,&an);
5643 }
5644 if(status<3){ /* if both ends are out of the frame do no further checks */
5645 if(status>0){
5646 p[nx]=p3;
5647 ix[nx]=i;
5648 jx[nx]=status;
5649 kx[nx]=2;
5650 ang[nx]=an;
5651 nx++;
5652 }
5653 status=CrossClipSide(&frameLt,angFrameLt,&p1,&p2,&p3,&an);
5654 }
5655 if(status<3){
5656 if(status>0){
5657 p[nx]=p3;
5658 ix[nx]=i;
5659 jx[nx]=status;
5660 kx[nx]=3;
5661 ang[nx]=an;
5662 nx++;
5663 }
5664 status=CrossClipSide(&frameBm,angFrameBm,&p1,&p2,&p3,&an);
5665 }
5666 if(status<3){
5667 if(status>0){
5668 p[nx]=p3;
5669 ix[nx]=i;
5670 jx[nx]=status;
5671 kx[nx]=4;
5672 ang[nx]=an;
5673 nx++;
5674 }
5675 }
5676 } /* for */
5677 if(nx==0){ /* no crossings */
5678 if(status==0) ; // the whole contour is inside the clip rect
5679 else{
5680 an=0.;
5681 for(i=0;i<j;i++){
5682 if(i==j-1) x=con[0].x-con[i].x;
5683 else x=con[i+1].x-con[i].x;
5684 if(x<-3.141592654) x=6.283185307+x;
5685 if(x>3.141592654) x=x-6.283185307;
5686 an+=x;
5687 }
5688 if(an<-6.){
5689 con[0]=RtBm;
5690 con[1]=TpRt;
5691 con[2]=LtTp;
5692 con[3]=BmLt;
5693 j=4;
5694 }
5695 else j=0;
5696 }
5697 }
5698 else{
5699 // if two crossings occur for a segment, 'enter' comes first and then comes 'leave'
5700 for(i=1;i<nx;i++){
5701 if(ix[i]==ix[i-1]){
5702 if(jx[i]==1){ // enter
5703 p3=p[i];
5704 p[i]=p[i-1];
5705 p[i-1]=p3;
5706 m=jx[i];
5707 jx[i]=jx[i-1];
5708 jx[i-1]=m;
5709 m=kx[i];
5710 kx[i]=kx[i-1];
5711 kx[i-1]=m;
5712 an=ang[i];
5713 ang[i]=ang[i-1];
5714 ang[i-1]=an;
5715 }
5716 }
5717 }
5718
5719 // add corner points after 'leave' and before 'enter'
5720 for(i=0;i<nx;i++){
5721 if(jx[i]==2){ // leave
5722 m=i+1;
5723 if(m==nx) m=0; // next enter point
5724 k=kx[i]; // side of leave point
5725 if(k!=kx[m]) do{
5726 k++;
5727 if(k==5) k=1;
5728 for(l=nx;l>i+1;l--){ // make room for a corner point
5729 p[l]=p[l-1];
5730 ix[l]=ix[l-1];
5731 jx[l]=jx[l-1];
5732 kx[l]=kx[l-1];
5733 ang[l]=ang[l-1];
5734 }
5735 switch(k){
5736 case 1:
5737 p[i+1]=RtBm;
5738 break;
5739 case 2:
5740 p[i+1]=TpRt;
5741 break;
5742 case 3:
5743 p[i+1]=LtTp;
5744 break;
5745 case 4:
5746 p[i+1]=BmLt;
5747 break;
5748 }
5749 ix[i+1]=ix[i];
5750 jx[i+1]=4; // status is corner point
5751 kx[i+1]=k;
5752 ang[i+1]=0.;
5753 nx++;
5754 i++;
5755 if(m!=0) m++;
5756 }while(k!=kx[m]); /* while */
5757 } /* if */
5758 } /* for */
5759 /* generate new contour */
5760 k=0; /* pending point */
5761 l=0; /* index into con1[] */
5762 m=0; /* index into p[] */
5763 n=0;
5764 for(i=0;i<nx;i++){ /* do for crossing points and corner points */
5765 if(jx[i]==2){ /* exit */
5766 for(o=k;o<=ix[i];o++) con1[l++]=con[o]; /* add preceeding points */
5767 con1[l++]=p[i];
5768 if(i==0) n=1;
5769 }
5770 else{ /* enter or corner */
5771 con1[l++]=p[i];
5772 k=ix[i]+1;
5773 }
5774 }
5775 if(n==1) for(o=k;o<j;o++) con1[l++]=con[o];
5776 for(i=0;i<l;i++) con[i]=con1[i];
5777 j=l;
5778 } /* else */
5779 } /* if */
5780 /* contour has been clipped, now we move on to line segments */
5781 k=blk->start;
5782 l=blk->pte;
5783 m=0;
5784 for(i=k;i<l;i++){
5785 p1.x=seg[i].sx;
5786 p1.y=1.570796327-seg[i].sy;
5787 p2.x=seg[i].ex;
5788 p2.y=1.570796327-seg[i].ey;
5789 nx=0;
5790 status=CrossClipSide(&frameRt,angFrameRt,&p1,&p2,&p3,&an);
5791 if(status<3){ // crossing or inside
5792 if(status>0){ // crossing
5793 p[nx]=p3; // crossing point in polar coordinate
5794 ix[nx]=i; // start point index
5795 jx[nx]=status; // status 1 or 2
5796 kx[nx]=1; // clip rect side - right
5797 ang[nx]=an; // distance from right-bottom corner
5798 nx++;
5799 }
5800 status=CrossClipSide(&frameTp,angFrameTp,&p1,&p2,&p3,&an);
5801 }
5802 if(status<3){ /* if both ends are out of the frame do no further checks */
5803 if(status>0){
5804 p[nx]=p3;
5805 ix[nx]=i;
5806 jx[nx]=status;
5807 kx[nx]=2;
5808 ang[nx]=an;
5809 nx++;
5810 }
5811 status=CrossClipSide(&frameLt,angFrameLt,&p1,&p2,&p3,&an);
5812 }
5813 if(status<3){
5814 if(status>0){
5815 p[nx]=p3;
5816 ix[nx]=i;
5817 jx[nx]=status;
5818 kx[nx]=3;
5819 ang[nx]=an;
5820 nx++;
5821 }
5822 status=CrossClipSide(&frameBm,angFrameBm,&p1,&p2,&p3,&an);
5823 }
5824 if(status<3){
5825 if(status>0){
5826 p[nx]=p3;
5827 ix[nx]=i;
5828 jx[nx]=status;
5829 kx[nx]=4;
5830 ang[nx]=an;
5831 nx++;
5832 }
5833 }
5834 if(nx==0){
5835 if(status==3); /* out of clipping rect */
5836 else seg1[m++]=seg[i];
5837 }
5838 else if(nx==1){
5839 if(jx[0]==1){ /* comming in */
5840 seg1[m].sx=p[0].x;
5841 seg1[m].sy=1.570796327-p[0].y;
5842 seg1[m].ex=seg[i].ex;
5843 seg1[m].ey=seg[i].ey;
5844 m++;
5845 }
5846 else{ /* going out */
5847 seg1[m].sx=seg[i].sx;
5848 seg1[m].sy=seg[i].sy;
5849 seg1[m].ex=p[0].x;
5850 seg1[m].ey=1.570796327-p[0].y;
5851 m++;
5852 }
5853 }
5854 else{
5855 seg1[m].sx=p[0].x;
5856 seg1[m].sy=1.570796327-p[0].y;
5857 seg1[m].ex=p[1].x;
5858 seg1[m].ey=1.570796327-p[1].y;
5859 m++;
5860 }
5861 }
5862 for(i=0;i<j;i++){
5863 seg[i].sx=con[i].x;
5864 seg[i].sy=1.570796327-con[i].y;
5865 if(i==j-1) k=0;
5866 else k=i+1;
5867 seg[i].ex=con[k].x;
5868 seg[i].ey=1.570796327-con[k].y;
5869 }
5870 cont->start=0;
5871 cont->pte=j;
5872 blk->start=j;
5873 for(i=0;i<m;i++) seg[j++]=seg1[i];
5874 blk->pte=j;
5875}
5876
5877/********************************************************************
5878 ToDB
5879
5880 generate stream data
5881********************************************************************/
5882char *CThreeD::ToDB(int *len)
5883{
5884 int i,j,k,l,m,n;
5885 CMPpart *p;
5886 Point2d *p2;
5887 Point3d *p3;
5888 char *e;
5889
5890 n=0;
5891 n+=4; /* npart */
5892 for(i=0;i<npart;i++){
5893 n+=4; /* i */
5894 p=part[i];
5895 n+=4; /* kind */
5896 j=strlen(p->desc.c_str());
5897 n+=4; /* j */
5898 n+=j+1; /* strcpy */
5899 n+=4; /* refPart */
5900 n+=8*6; /* x0,y0,z0,theta0,phi0,psi0 */
5901 n+=4; /* color */
5902 n+=4; /* showDrawing */
5903 switch(p->kind){
5904 case MPpoint:
5905 break;
5906 case MPlines:
5907 k=reinterpret_cast<CMPlines *>(p)->npoint;
5908 n+=4; /* k */
5909 n+=k*24; /* x,y,z */
5910 n+=4; /* line */
5911 break;
5912 case MPplane:
5913 k=reinterpret_cast<CMPplane *>(p)->npoint;
5914 n+=4; /* k */
5915 n+=k*16; /* x,y */
5916 n+=4; /* isPolyhed */
5917 n+=4; /* isGroup */
5918 break;
5919 case MPbox:
5920 n+=24; /* lx,ly,lz */
5921 n+=4; /* visibleSurface */
5922 break;
5923 case MPsphere:
5924 n+=8; /* r */
5925 break;
5926 case MPbodyOfRot:
5927 n+=4; /* npoint */
5928 n+=32; /* p0.x,p0.y,p1.x,p1.y */
5929 n+=16; /* isTube,oblong,appear,visibleSurfac */
5930 break;
5931 }
5932 }
5933 *len=n;
5934 e=(char *)malloc(sizeof(char)*n);
5935 n=0;
5936 *(int *)(e+n)=npart;
5937 n+=4;
5938 for(i=0;i<npart;i++){
5939 *(int *)(e+n)=i;
5940 n+=4;
5941 p=part[i];
5942 *(int *)(e+n)=p->kind;
5943 n+=4;
5944 j=strlen(p->desc.c_str());
5945 *(int *)(e+n)=j;
5946 n+=4;
5947 strcpy(e+n,p->desc.c_str());
5948 n+=j+1;
5949 *(int *)(e+n)=p->refPart;
5950 n+=4;
5951 *(double *)(e+n)=p->x0;
5952 n+=8;
5953 *(double *)(e+n)=p->y0;
5954 n+=8;
5955 *(double *)(e+n)=p->z0;
5956 n+=8;
5957 *(double *)(e+n)=p->theta0;
5958 n+=8;
5959 *(double *)(e+n)=p->phi0;
5960 n+=8;
5961 *(double *)(e+n)=p->psi0;
5962 n+=8;
5963 *(int *)(e+n)=p->color;
5964 n+=4;
5965 *(int *)(e+n)=p->showDrawing;
5966 n+=4;
5967 switch(p->kind){
5968 case MPpoint:
5969 break;
5970 case MPlines:
5971 k=reinterpret_cast<CMPlines *>(p)->npoint;
5972 *(int *)(e+n)=k;
5973 n+=4;
5974 for(j=0;j<k;j++){
5975 p3=&reinterpret_cast<CMPlines *>(p)->pt[j];
5976 *(double *)(e+n)=p3->x;
5977 n+=8;
5978 *(double *)(e+n)=p3->y;
5979 n+=8;
5980 *(double *)(e+n)=p3->z;
5981 n+=8;
5982 }
5983 *(int *)(e+n)=reinterpret_cast<CMPlines *>(p)->line;
5984 n+=4;
5985 break;
5986 case MPplane:
5987 k=reinterpret_cast<CMPplane *>(p)->npoint;
5988 *(int *)(e+n)=k;
5989 n+=4;
5990 for(j=0;j<k;j++){
5991 p2=&reinterpret_cast<CMPplane *>(p)->pt[j];
5992 *(double *)(e+n)=p2->x;
5993 n+=8;
5994 *(double *)(e+n)=p2->y;
5995 n+=8;
5996 }
5997 *(int *)(e+n)=reinterpret_cast<CMPplane *>(p)->isPolyhed;
5998 n+=4;
5999 *(int *)(e+n)=reinterpret_cast<CMPplane *>(p)->isGroup;
6000 n+=4;
6001 break;
6002 case MPbox:
6003 *(double *)(e+n)=reinterpret_cast<CMPbox *>(p)->lx;
6004 n+=8;
6005 *(double *)(e+n)=reinterpret_cast<CMPbox *>(p)->ly;
6006 n+=8;
6007 *(double *)(e+n)=reinterpret_cast<CMPbox *>(p)->lz;
6008 n+=8;
6009 *(int *)(e+n)=reinterpret_cast<CMPbox *>(p)->visibleSurface;
6010 n+=4;
6011 break;
6012 case MPsphere:
6013 *(double *)(e+n)=reinterpret_cast<CMPsphere *>(p)->r;
6014 n+=8;
6015 break;
6016 case MPbodyOfRot:
6017 *(int *)(e+n)=reinterpret_cast<CMPbodyOfRot *>(p)->npoint;
6018 n+=4;
6019 *(double *)(e+n)=reinterpret_cast<CMPbodyOfRot *>(p)->p0.x;
6020 n+=8;
6021 *(double *)(e+n)=reinterpret_cast<CMPbodyOfRot *>(p)->p0.y;
6022 n+=8;
6023 *(double *)(e+n)=reinterpret_cast<CMPbodyOfRot *>(p)->p1.x;
6024 n+=8;
6025 *(double *)(e+n)=reinterpret_cast<CMPbodyOfRot *>(p)->p1.y;
6026 n+=8;
6027 *(int *)(e+n)=reinterpret_cast<CMPbodyOfRot *>(p)->isTube;
6028 n+=4;
6029 *(int *)(e+n)=reinterpret_cast<CMPbodyOfRot *>(p)->oblong;
6030 n+=4;
6031 *(int *)(e+n)=reinterpret_cast<CMPbodyOfRot *>(p)->appear;
6032 n+=4;
6033 *(int *)(e+n)=reinterpret_cast<CMPbodyOfRot *>(p)->visibleSurface;
6034 n+=4;
6035 break;
6036 }
6037 }
6038 return e;
6039}
6040
6041/********************************************************************
6042 FromDB
6043
6044 load stream data
6045********************************************************************/
6046void CThreeD::FromDB(int len,char *e,double loc[6])
6047{
6048 int i,j,k,l,m,n,i1,i2;
6049 char s[200];
6050 CMPpart *p;
6051 Point2d *p2;
6052 Point3d *p3;
6053
6054 n=0;
6055 l=npart;
6056 m=*(int *)e;
6057 n+=4;
6058 for(i=0;i<m;i++){
6059 /* *(int *)(e+n)=i */
6060 n+=4;
6061 i1=*(int *)(e+n); // kind
6062 n+=4;
6063 i2=*(int *)(e+n); // strlen
6064 n+=4;
6065 strcpy(s,e+n);
6066 n+=i2+1;
6067 i2=*(int *)(e+n)+l; // refPart
6068 if(i==0) i2=0;
6069 n+=4;
6070 switch(i1){
6071 case MPpoint:
6072 part.push_back(new CMPpoint(s,i2,i1));
6073 break;
6074 case MPlines:
6075 part.push_back(new CMPlines(s,i2,i1));
6076 break;
6077 case MPplane:
6078 part.push_back(new CMPplane(s,i2,i1));
6079 break;
6080 case MPbox:
6081 part.push_back(new CMPbox(s,i2,i1));
6082 break;
6083 case MPsphere:
6084 part.push_back(new CMPsphere(s,i2,i1));
6085 break;
6086 case MPbodyOfRot:
6087 part.push_back(new CMPbodyOfRot(s,i2,i1));
6088 break;
6089 default:
6090 break;
6091 }
6092 npart++;
6093 p=part[npart-1];
6094 if(i==0){
6095 p->x0=loc[0];
6096 p->y0=loc[1];
6097 p->z0=loc[2];
6098 p->theta0=loc[3];
6099 p->phi0=loc[4];
6100 p->psi0=loc[5];
6101 n+=48;
6102 }
6103 else{
6104 p->x0=*(double *)(e+n);
6105 n+=8;
6106 p->y0=*(double *)(e+n);
6107 n+=8;
6108 p->z0=*(double *)(e+n);
6109 n+=8;
6110 p->theta0=*(double *)(e+n);
6111 n+=8;
6112 p->phi0=*(double *)(e+n);
6113 n+=8;
6114 p->psi0=*(double *)(e+n);
6115 n+=8;
6116 }
6117 p->color=*(int *)(e+n);
6118 n+=4;
6119 p->showDrawing=*(int *)(e+n);
6120 n+=4;
6121 switch(p->kind){
6122 case MPpoint:
6123 break;
6124 case MPlines:
6125 k=*(int *)(e+n);
6126 n+=4;
6127 for(j=0;j<k;j++){
6128 reinterpret_cast<CMPlines *>(p)->SetPoint(*(double *)(e+n),*(double *)(e+n+8),*(double *)(e+n+16));
6129 n+=24;
6130 }
6131 reinterpret_cast<CMPlines *>(p)->line=*(int *)(e+n);
6132 n+=4;
6133 break;
6134 case MPplane:
6135 k=*(int *)(e+n);
6136 n+=4;
6137 for(j=0;j<k;j++){
6138 reinterpret_cast<CMPplane *>(p)->SetPoint(*(double *)(e+n),*(double *)(e+n+8));
6139 n+=16;
6140 }
6141 reinterpret_cast<CMPplane *>(p)->isPolyhed=*(int *)(e+n);
6142 n+=4;
6143 reinterpret_cast<CMPplane *>(p)->isGroup=*(int *)(e+n);
6144 n+=4;
6145 break;
6146 case MPbox:
6147 reinterpret_cast<CMPbox *>(p)->SetSize(*(double *)(e+n),*(double *)(e+n+8),*(double *)(e+n+16));
6148 n+=24;
6149 reinterpret_cast<CMPbox *>(p)->visibleSurface=*(int *)(e+n);
6150 n+=4;
6151 break;
6152 case MPsphere:
6153 reinterpret_cast<CMPsphere *>(p)->r=*(double *)(e+n);
6154 n+=8;
6155 break;
6156 case MPbodyOfRot:
6157 reinterpret_cast<CMPbodyOfRot *>(p)->npoint=*(int *)(e+n);
6158 n+=4;
6159 reinterpret_cast<CMPbodyOfRot *>(p)->p0.x=*(double *)(e+n);
6160 n+=8;
6161 reinterpret_cast<CMPbodyOfRot *>(p)->p0.y=*(double *)(e+n);
6162 n+=8;
6163 reinterpret_cast<CMPbodyOfRot *>(p)->p1.x=*(double *)(e+n);
6164 n+=8;
6165 reinterpret_cast<CMPbodyOfRot *>(p)->p1.y=*(double *)(e+n);
6166 n+=8;
6167 reinterpret_cast<CMPbodyOfRot *>(p)->isTube=*(int *)(e+n);
6168 n+=4;
6169 reinterpret_cast<CMPbodyOfRot *>(p)->oblong=*(int *)(e+n);
6170 n+=4;
6171 reinterpret_cast<CMPbodyOfRot *>(p)->appear=*(int *)(e+n);
6172 n+=4;
6173 reinterpret_cast<CMPbodyOfRot *>(p)->visibleSurface=*(int *)(e+n);
6174 n+=4;
6175 break;
6176 }
6177 }
6178}
6179
6180/******************************************************************************
6181 CFarNear Far-Near relationship of 3D Drawing Objects
6182******************************************************************************/
6183
6184
6185CFarNear::CFarNear(int numpart1)
6186 : npart1(numpart1)
6187{
6188 int i;
6189
6190 nmax=npart1*(npart1-1)/2;
6191 farNear=new StartPTE[nmax+1]; // (StartPTE *)NewPtr(sizeof(StartPTE)*nmax);
6192 farNearIndex=new long[npart1+1]; // (long *)NewPtr(sizeof(long)*(npart1+1));
6193 for(i=0;i<=npart1;i++) farNearIndex[i]=0L;
6194 n=0;
6195}
6196
6197/******************************************************************************
6198 Dispose
6199
6200******************************************************************************/
6201
6202CFarNear::~CFarNear()
6203{
6204 delete[] farNear; // DisposPtr((Ptr)farNear);
6205 delete[] farNearIndex; // DisposPtr((Ptr)farNearIndex);
6206}
6207
6208/******************************************************************************
6209 AddFarNear
6210
6211******************************************************************************/
6212
6213void CFarNear::AddFarNear(int far,int near)
6214{
6215 StartPTE a;
6216 long m;
6217 int i;
6218
6219 a.start=far;
6220 a.pte=near;
6221 m=farNearIndex[a.start+1];
6222 for(i=a.start+1;i<=npart1;i++) farNearIndex[i]++;
6223 // BlockMove((Ptr)&farNear[m],(Ptr)&farNear[m+1],(n-m)*sizeof(StartPTE));
6224 for(i=n;i>=m;i--) farNear[i+1]=farNear[i];
6225 farNear[m]=a;
6226 n++;
6227 // outstr("in AddFarNear()\n");
6228 // for(i=0;i<n;i++){
6229 // sprintf(msg," i=%i farNear[i].start=%i farNear[i].pte=%i \n",i,farNear[i].start,farNear[i].pte);
6230 // outstr(msg);
6231 // }
6232 // flush();
6233
6234}
6235
6236/******************************************************************************
6237 inferFarNear
6238
6239******************************************************************************/
6240
6241bool CFarNear::InferFarNear(int p1,int p2,int *near)
6242{
6243 // test if p1 is farther
6244 *near=p2;
6245 if(RecurseFarNear(p1,p2,0)) return true;
6246 // test if p2 is farther
6247 *near=p1;
6248 if(RecurseFarNear(p2,p1,0)) return true;
6249 return false;
6250}
6251
6252/******************************************************************************
6253 recurseFarNear
6254
6255******************************************************************************/
6256
6257bool CFarNear::RecurseFarNear(int p1,int p2,long depth)
6258{
6259 long i,j,k;
6260
6261 depth++;
6262 if(depth>n) return false; // escape loop
6263
6264 j=farNearIndex[p1];
6265 k=farNearIndex[p1+1];
6266 for(i=j;i<k;i++){
6267 if(farNear[i].pte==p2) return true;
6268 if(RecurseFarNear(farNear[i].pte,p2,depth)) return true;
6269 }
6270 return false;
6271}
6272
6273/******************************************************************************
6274 SetSorted
6275
6276******************************************************************************/
6277
6278int CFarNear::SetSorted(int *sorted)
6279{
6280 int i,j,l;
6281 long k,m;
6282 int *status;
6283
6284 // outstr("in SetSorted()\n");
6285 // sprintf(msg," npart1=%i n=%i \n",npart1,n);
6286 // outstr(msg);
6287 // flush();
6288 // for(i=0;i<n;i++){
6289 // sprintf(msg," farNear[i].start=%5i farNear[i].pte=%5i\n",farNear[i].start,farNear[i].pte);
6290 // outstr(msg);
6291 // }
6292 // flush();
6293
6294 status=new int[npart1]; // (int *)NewPtr(sizeof(short)*npart1);
6295 for(i=0;i<npart1;i++) status[i]=0;
6296 j=0;
6297 i=0;
6298 do{
6299 for(k=0;k<n;k++){
6300 l=farNear[k].start;
6301 if(l<10000){
6302 switch (status[l]){
6303 case 0:
6304 for(m=0;m<n;m++) if(l==farNear[m].pte) break;
6305 if(m==n){
6306 sorted[j++]=l;
6307 status[l]=2;
6308 farNear[k].start=10000;
6309 farNear[k].pte=10000;
6310 }
6311 else status[l]=1;
6312 break;
6313 case 1:
6314 break;
6315 case 2:
6316 case 3:
6317 farNear[k].start=10000;
6318 farNear[k].pte=10000;
6319 break;
6320 default:
6321 break;
6322 }
6323 }
6324 }
6325 m=0;
6326 for(k=0;k<npart1;k++){
6327 l=status[k];
6328 if(l==1){
6329 m++;
6330 status[k]=0;
6331 }
6332 if(l==2){
6333 if(m<10000) m+=10000;
6334 status[k]=3;
6335 }
6336 }
6337 if(m<10000) i++;
6338 } while(i<2);
6339 for(k=0;k<npart1;k++) if(status[k]!=3) sorted[j++]=k;
6340 delete[] status; // DisposPtr((Ptr)status);
6341 return m;
6342}
6343
6344/**************************************************
6345
6346 Initiate adding a rigid body element
6347
6348 should be followed by following calls as required
6349 setlocal
6350 setweight
6351 setdensity
6352 setnpoint
6353 set1Dpoint
6354 set2Dpoint
6355 set3Dpoint
6356
6357***************************************************/
6358
6359int addelm(char *s,int ref,int knd1)
6360{
6361 int knd,piggy;
6362 piggy=knd1/100;
6363 knd=knd1-piggy*100;
6364 switch(knd){
6365 case MPpoint:
6366 thd.part.push_back(new CMPpoint(s,ref,knd));
6367 break;
6368 case MPlines:
6369 thd.part.push_back(new CMPlines(s,ref,knd));
6370 break;
6371 case MPplane:
6372 thd.part.push_back(new CMPplane(s,ref,knd));
6373 if(piggy) (reinterpret_cast<CMPplane *>(thd.part[thd.npart]))->isGroup=1;
6374 break;
6375 case MPbox:
6376 thd.part.push_back(new CMPbox(s,ref,knd));
6377 if(piggy) (reinterpret_cast<CMPbox *>(thd.part[thd.npart]))->objid=piggy;
6378 break;
6379 case MPsphere:
6380 thd.part.push_back(new CMPsphere(s,ref,knd));
6381 break;
6382 case MPbodyOfRot:
6383 thd.part.push_back(new CMPbodyOfRot(s,ref,knd));
6384 break;
6385 default:
6386 return -1;
6387 break;
6388 }
6389 thd.npart++;
6390 return thd.npart-1;
6391}
6392
6393/****************************************************
6394
6395 set local coordinate
6396
6397****************************************************/
6398void setlocal(double x,double y,double z,double theta,double phi,double psi)
6399{
6400 thd.part[thd.npart-1]->x0=x;
6401 thd.part[thd.npart-1]->y0=y;
6402 thd.part[thd.npart-1]->z0=z;
6403 thd.part[thd.npart-1]->theta0=theta;
6404 thd.part[thd.npart-1]->phi0=phi;
6405 thd.part[thd.npart-1]->psi0=psi;
6406}
6407
6408/***************************************************
6409
6410 set element color
6411
6412***************************************************/
6413void setcolor(int c)
6414{
6415 thd.part[thd.npart-1]->color=c;
6416}
6417
6418/***************************************************
6419
6420 set line appearance
6421
6422***************************************************/
6423int setline(int c)
6424{
6425 if(thd.part[thd.npart-1]->kind!=MPlines) return -1;
6426 (reinterpret_cast<CMPlines *>(thd.part[thd.npart-1]))->line=c;
6427 return 0;
6428}
6429
6430/***************************************************
6431
6432 set if tube for a body of rotation
6433
6434***************************************************/
6435int isTube()
6436{
6437 if(thd.part[thd.npart-1]->kind!=MPbodyOfRot) return -1;
6438 (reinterpret_cast<CMPbodyOfRot *>(thd.part[thd.npart-1]))->isTube=true;
6439 return 0;
6440}
6441
6442/***************************************************
6443
6444 set one dimensional point as in sphere radius
6445
6446***************************************************/
6447int set1Dpoint(double x)
6448{
6449 switch(thd.part[thd.npart-1]->kind){
6450 case MPsphere:
6451 (reinterpret_cast<CMPsphere *>(thd.part[thd.npart-1]))->SetRadius(x);
6452 break;
6453 default:
6454 return -1;
6455 break;
6456 }
6457 return 0;
6458}
6459
6460/***************************************************
6461
6462 set two dimensional point as in an apex of plane
6463
6464***************************************************/
6465int set2Dpoint(double x,double y)
6466{
6467 switch(thd.part[thd.npart-1]->kind){
6468 case MPplane:
6469 (reinterpret_cast<CMPplane *>(thd.part[thd.npart-1]))->SetPoint(x,y);
6470 break;
6471 case MPbodyOfRot:
6472 (reinterpret_cast<CMPbodyOfRot *>(thd.part[thd.npart-1]))->SetPoint(x,y);
6473 break;
6474 default:
6475 return -1;
6476 break;
6477 }
6478 return 0;
6479}
6480
6481/***************************************************
6482
6483 set three dimensional point as in a point in line
6484
6485***************************************************/
6486int set3Dpoint(double x,double y,double z)
6487{
6488 switch(thd.part[thd.npart-1]->kind){
6489 case MPlines:
6490 (reinterpret_cast<CMPlines *>(thd.part[thd.npart-1]))->SetPoint(x,y,z);
6491 break;
6492 case MPbox:
6493 (reinterpret_cast<CMPbox *>(thd.part[thd.npart-1]))->SetSize(x,y,z);
6494 break;
6495 default:
6496 return -1;
6497 break;
6498 }
6499 return 0;
6500}
6501
6502void viewFrom(double a,double b,double r)
6503{
6504 thd.Azi=a;
6505 thd.Ele=b;
6506 thd.Dst=r;
6507}
6508
6509void drawMessage(int i)
6510{
6511 thd.message=i;
6512}
6513
6514void postscript3D(int i)
6515{
6516 thd.postscript=i;
6517}
6518
6519void viewFixed(double tgx,double tgy,double tgz,int foc)
6520{
6521 thd.contAll=false;
6522 thd.Tgx=tgx;
6523 thd.Tgy=tgy;
6524 thd.Tgz=tgz;
6525 thd.Foc=foc/2;
6526}
6527
6528void omnidirect(int omni)
6529{
6530 thd.omni=omni;
6531}
6532
6533int draw3D(int wid,int hei)
6534{
6535 thd.Wid=wid/2;
6536 thd.Hei=hei/2;
6537 thd.Draw();
6538 return 0;
6539}
6540
6541void getCoord(int ref,double x,double y,double z,int *ix,int *iy)
6542{
6543 /* to be implemented some time later */
6544}
6545
6546char *insert3D(int *len)
6547{
6548 char *s;
6549 s=thd.ToDB(len);
6550}
6551
6552void select3D(int len,char *d,double x,double y,double z,double theta,double phi,double psi)
6553{
6554 double loc[6];
6555 loc[0]=x;
6556 loc[1]=y;
6557 loc[2]=z;
6558 loc[3]=theta;
6559 loc[4]=phi;
6560 loc[5]=psi;
6561 thd.FromDB(len,d,loc);
6562}
6563
6564void parts3D(int style)
6565{
6566 char s[200];
6567 int i,j,k,l;
6568 CMPpart *p;
6569 Point2d *p2;
6570 Point3d *p3;
6571 if(style==1){
6572 for(i=1;i<thd.npart;i++){
6573 sprintf(s," addelm(\"plane %i\",0,2);\n",i-1);
6574 outstr(s);
6575 p=thd.part[i];
6576 sprintf(s," setlocal(%g,%g,%g,%g,%g,%g);\n",p->x0,p->y0,p->z0,p->theta0,p->phi0,p->psi0);
6577 outstr(s);
6578 k=reinterpret_cast<CMPplane *>(p)->npoint;
6579 for(j=0;j<k;j++){
6580 p2=&reinterpret_cast<CMPplane *>(p)->pt[j];
6581 sprintf(s," set2Dpoint(%g,%g);\n",p2->x,p2->y);
6582 outstr(s);
6583 }
6584 outstr(" setcolor(5);\n");
6585 flush();
6586 }
6587 return;
6588 }
6589 sprintf(s,"thd.npart=%i\n",thd.npart);
6590 outstr(s);
6591 for(i=0;i<thd.npart;i++){
6592 sprintf(s,"part index=%i ",i);
6593 outstr(s);
6594 p=thd.part[i];
6595 sprintf(s,"kind=%i ",p->kind);
6596 outstr(s);
6597 strcpy(s,p->desc.c_str());
6598 outstr(s);
6599 sprintf(s,"\n refPart=%i ",p->refPart);
6600 outstr(s);
6601 sprintf(s,"x0=%g ",p->x0);
6602 outstr(s);
6603 sprintf(s,"y0=%g ",p->y0);
6604 outstr(s);
6605 sprintf(s,"z0=%g ",p->z0);
6606 outstr(s);
6607 sprintf(s,"theta0=%g ",p->theta0);
6608 outstr(s);
6609 sprintf(s,"phi0=%g ",p->phi0);
6610 outstr(s);
6611 sprintf(s,"psi0=%g ",p->psi0);
6612 outstr(s);
6613 sprintf(s,"color=%i ",p->color);
6614 outstr(s);
6615 sprintf(s,"showDrawing=%i\n",p->showDrawing);
6616 outstr(s);
6617 switch(p->kind){
6618 case MPpoint:
6619 break;
6620 case MPlines:
6621 k=reinterpret_cast<CMPlines *>(p)->npoint;
6622 sprintf(s,"npoint=%i\n",k);
6623 for(j=0;j<k;j++){
6624 p3=&reinterpret_cast<CMPlines *>(p)->pt[j];
6625 sprintf(s,"x=%g y=%g z=%g\n",p3->x,p3->y,p3->z);
6626 outstr(s);
6627 }
6628 sprintf(s,"line=%i\n",reinterpret_cast<CMPlines *>(p)->line);
6629 outstr(s);
6630 break;
6631 case MPplane:
6632 k=reinterpret_cast<CMPplane *>(p)->npoint;
6633 sprintf(s,"npoint=%i\n",k);
6634 for(j=0;j<k;j++){
6635 p2=&reinterpret_cast<CMPplane *>(p)->pt[j];
6636 sprintf(s,"x=%g y=%g\n",p2->x,p2->y);
6637 outstr(s);
6638 }
6639 sprintf(s,"isPolyhed=%i isGroup=%i\n",reinterpret_cast<CMPplane *>(p)->isPolyhed,
6640 reinterpret_cast<CMPplane *>(p)->isGroup);
6641 break;
6642 case MPbox:
6643 sprintf(s,"lx=%g ly=%g lz=%g visibleSurface=%i\n",
6644 reinterpret_cast<CMPbox *>(p)->lx,reinterpret_cast<CMPbox *>(p)->ly,
6645 reinterpret_cast<CMPbox *>(p)->lz,
6646 reinterpret_cast<CMPbox *>(p)->visibleSurface);
6647 outstr(s);
6648 break;
6649 case MPsphere:
6650 sprintf(s,"r=%g\n",reinterpret_cast<CMPsphere *>(p)->r);
6651 outstr(s);
6652 break;
6653 case MPbodyOfRot:
6654 sprintf(s,"npoint=%i ",reinterpret_cast<CMPbodyOfRot *>(p)->npoint);
6655 outstr(s);
6656 sprintf(s,"p0.x=%g p0.y=%g ",reinterpret_cast<CMPbodyOfRot *>(p)->p0.x,
6657 reinterpret_cast<CMPbodyOfRot *>(p)->p0.y);
6658 outstr(s);
6659 sprintf(s,"p1.x=%g p1.y=%g\n",reinterpret_cast<CMPbodyOfRot *>(p)->p1.x,
6660 reinterpret_cast<CMPbodyOfRot *>(p)->p1.y);
6661 outstr(s);
6662 sprintf(s,"isTube=%i ",reinterpret_cast<CMPbodyOfRot *>(p)->isTube);
6663 outstr(s);
6664 sprintf(s,"oblong=%i ",reinterpret_cast<CMPbodyOfRot *>(p)->oblong);
6665 outstr(s);
6666 sprintf(s,"appear=%i ",reinterpret_cast<CMPbodyOfRot *>(p)->appear);
6667 outstr(s);
6668 sprintf(s,"visibleSurface=%i\n",reinterpret_cast<CMPbodyOfRot *>(p)->visibleSurface);
6669 outstr(s);
6670 }
6671 }
6672 flush();
6673}
6674