/******************************************************************** * Description: * Author: George-Catalin Serbanut * * Copyright (c) 2005 George-Catalin Serbanut All rights reserved. * ********************************************************************/ /* * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ #include #include #include #include #include #include "include/data_igs.hh" #include "include/converting.hh" #include "include/math.hh" #include "include/algebra.hh" #include "include/trigonometry.hh" #include "include/point.hh" #include "include/line.hh" #include "include/plan.hh" #include "include/crystal.hh" using namespace std; Convert::Convert() { } Convert::~Convert() { } void Convert::multiplication() { vector ldigs; DataIGS ldigs1; double ra = 22.5*acos(-1.0)/180.0; double x1, y1, z1; double x2, y2, z2; for(int i=0; i<16; i++) for(int j=0; j(i))+y1*sin(ra*static_cast(i)); y2 = y1*cos(ra*static_cast(i))-x1*sin(ra*static_cast(i)); z2 = z1; ldigs1.x[k] = x2; ldigs1.y[k] = y2; ldigs1.z[k] = z2; }; ldigs1.module = Convert::igsdata[j].module; ldigs1.row = Convert::igsdata[j].row; ldigs1.crystal = 10*i+Convert::igsdata[j].crystal; ldigs.push_back(ldigs1); }; Convert::igsdata.clear(); for(int i=0; i plan1, plan2; vector median; Point p; Crystal crys; fstream file; double deviation; for(int i=0; i0.0)||(math.Precision(median[2].GetZ()-1.0,math.precision)>0.0)) // cout<<"Warning: : "<p1) // if((p3!=p2)&&(p3!=p1)) // if((p4>p3)&&(p4!=p2)&&(p4!=p1)) { // double angle_line = 180.0*trig.Angle(Line(plan2[p1],plan2[p2]),Line(plan2[p3],plan2[p4]))/acos(-1.0); // if((math.Precision(10.0-math.Abs(angle_line),6)>=0.0)||(math.Precision(math.Abs(angle_line)-170.0,6)>=0.0)) // cout<<"Angle(("<(acos(-1.0)/2.0)) { Convert::perp = crys.GetPerp()*(-1.0); } else { Convert::perp = crys.GetPerp(); }; Convert::substract(); Convert::rotate_points(); Convert::pDz(); Convert::pTheta(); Convert::pPhi(); for(int j=0; j<8; j++) Convert::pts[j] = -1; Convert::compute_plane("lower"); Convert::compute_plane("upper"); // just numerical imprecision Convert::dg4.pAlp1 = (Convert::dg4.pAlp1+Convert::dg4.pAlp2)/2.0; if(Convert::math.Abs(Convert::dg4.pAlp1)<(1e-5)) Convert::dg4.pAlp1 = 0.0; Convert::dg4.pAlp2 = Convert::dg4.pAlp1; if(Convert::dg4.module==2) { Convert::dg4.pAlp1 *= -1.0; Convert::dg4.pAlp2 *= -1.0; }; Convert::g4data.push_back(Convert::dg4); }; return; } // rotating the axes - no need! void Convert::rotate_axes(int igs_contor) { for(int j=0; j<8; j++) { // Convert::digs.x[j] = Convert::igsdata[igs_contor].z[j]; // Convert::digs.y[j] = Convert::igsdata[igs_contor].x[j]; // Convert::digs.z[j] = Convert::igsdata[igs_contor].y[j]; Convert::digs.x[j] = Convert::igsdata[igs_contor].x[j]; Convert::digs.y[j] = Convert::igsdata[igs_contor].y[j]; Convert::digs.z[j] = Convert::igsdata[igs_contor].z[j]; }; return; } // ordering over plans: p1 is the lower one (the first 4 points) and // p2 is the upper one (the last 4 points) // in this version it can be found in the convert_igs_g4 function void Convert::order() { return; } // computing the center of the crystal position void Convert::center_position() { Convert::xc1 = 0.0; Convert::yc1 = 0.0; Convert::zc1 = 0.0; Convert::xc2 = 0.0; Convert::yc2 = 0.0; Convert::zc2 = 0.0; for(int j=0; j<4; j++) { Convert::xc1 += Convert::digs.x[j]; Convert::yc1 += Convert::digs.y[j]; Convert::zc1 += Convert::digs.z[j]; }; for(int j=4; j<8; j++) { Convert::xc2 += Convert::digs.x[j]; Convert::yc2 += Convert::digs.y[j]; Convert::zc2 += Convert::digs.z[j]; }; Convert::xc1 = Convert::xc1/4.0; Convert::yc1 = Convert::yc1/4.0; Convert::zc1 = Convert::zc1/4.0; Convert::xc2 = Convert::xc2/4.0; Convert::yc2 = Convert::yc2/4.0; Convert::zc2 = Convert::zc2/4.0; Convert::dg4.posX = (Convert::xc1+Convert::xc2)/2.0; Convert::dg4.posY = (Convert::yc1+Convert::yc2)/2.0; Convert::dg4.posZ = (Convert::zc1+Convert::zc2)/2.0; Convert::xc = Convert::dg4.posX; Convert::yc = Convert::dg4.posY; Convert::zc = Convert::dg4.posZ; return; } void Convert::substract() { for(int i=0; i<8; i++) { Convert::digs.x[i] -= Convert::dg4.posX; Convert::digs.y[i] -= Convert::dg4.posY; Convert::digs.z[i] -= Convert::dg4.posZ; }; Convert::xc1 -= Convert::dg4.posX; Convert::xc2 -= Convert::dg4.posX; Convert::yc1 -= Convert::dg4.posY; Convert::yc2 -= Convert::dg4.posY; Convert::zc1 -= Convert::dg4.posZ; Convert::zc2 -= Convert::dg4.posZ; return; } void Convert::rotate_points() { double x1, y1, z1; // before double x2, y2, z2; // after double theta, phi, tau; Angles angs; Point pt; Algebra alg; angs = Angles(Convert::perp.GetPoint()); theta = angs.GetTheta(); phi = angs.GetOmega(); // cout<<"Angles: "<Convert::digs.y[3]) { for(int i=0; i<8; i++) { x1 = Convert::digs.x[i]; y1 = Convert::digs.y[i]; z1 = Convert::digs.z[i]; // Oz ccw rotation x2 = -1.0*x1; y2 = -1.0*y1; z2 = z1; Convert::digs.x[i] = x2; Convert::digs.y[i] = y2; Convert::digs.z[i] = z2; }; // rotating planes centers for(int i=0; i<2; i++) { if(i==0) { x1 = Convert::xc1; y1 = Convert::yc1; z1 = Convert::zc1; } else { x1 = Convert::xc2; y1 = Convert::yc2; z1 = Convert::zc2; }; // Oz ccw rotation x2 = -1.0*x1; y2 = -1.0*y1; z2 = z1; if(i==0) { Convert::xc1 = x2; Convert::yc1 = y2; Convert::zc1 = z2; } else { Convert::xc2 = x2; Convert::yc2 = y2; Convert::zc2 = z2; }; }; tau -= acos(-1.0); }; Convert::dg4.theta = 180.0*theta/acos(-1.0); Convert::dg4.phi = 180.0*phi/acos(-1.0); Convert::dg4.tau = 180.0*tau/acos(-1.0); //cout<pDy) { Convert::pts[min] = l1; Convert::pts[min+1] = l2; pDy = math.Abs(yc-(Convert::digs.y[l1]+Convert::digs.y[l2])/2.0); }; Convert::pts[min+2] = -1; for(int j=min; j