Changeset 643e76 for src/vector_ops.cpp
- Timestamp:
- May 27, 2010, 2:24:34 PM (16 years ago)
- Branches:
- Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, Candidate_v1.7.0, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
- Children:
- 42a101
- Parents:
- 45ef76
- File:
-
- 1 edited
-
src/vector_ops.cpp (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
src/vector_ops.cpp
r45ef76 r643e76 143 143 return res; 144 144 }; 145 146 /** Calculates the intersection of the two lines that are both on the same plane.147 * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html148 * \param *out output stream for debugging149 * \param *Line1a first vector of first line150 * \param *Line1b second vector of first line151 * \param *Line2a first vector of second line152 * \param *Line2b second vector of second line153 * \return true - \a this will contain the intersection on return, false - lines are parallel154 */155 Vector GetIntersectionOfTwoLinesOnPlane(const Vector &Line1a, const Vector &Line1b, const Vector &Line2a, const Vector &Line2b)156 {157 Info FunctionInfo(__func__);158 159 Vector res;160 161 auto_ptr<GSLMatrix> M = auto_ptr<GSLMatrix>(new GSLMatrix(4,4));162 163 M->SetAll(1.);164 for (int i=0;i<3;i++) {165 M->Set(0, i, Line1a[i]);166 M->Set(1, i, Line1b[i]);167 M->Set(2, i, Line2a[i]);168 M->Set(3, i, Line2b[i]);169 }170 171 //Log() << Verbose(1) << "Coefficent matrix is:" << endl;172 //for (int i=0;i<4;i++) {173 // for (int j=0;j<4;j++)174 // cout << "\t" << M->Get(i,j);175 // cout << endl;176 //}177 if (fabs(M->Determinant()) > MYEPSILON) {178 Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;179 throw SkewException(__FILE__,__LINE__);180 }181 182 Log() << Verbose(1) << "INFO: Line1a = " << Line1a << ", Line1b = " << Line1b << ", Line2a = " << Line2a << ", Line2b = " << Line2b << "." << endl;183 184 185 // constuct a,b,c186 Vector a = Line1b - Line1a;187 Vector b = Line2b - Line2a;188 Vector c = Line2a - Line1a;189 Vector d = Line2b - Line1b;190 Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;191 if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {192 res.Zero();193 Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;194 throw LinearDependenceException(__FILE__,__LINE__);195 }196 197 // check for parallelity198 Vector parallel;199 double factor = 0.;200 if (fabs(a.ScalarProduct(b)*a.ScalarProduct(b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {201 parallel = Line1a - Line2a;202 factor = parallel.ScalarProduct(a)/a.Norm();203 if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {204 res = Line2a;205 Log() << Verbose(1) << "Lines conincide." << endl;206 return res;207 } else {208 parallel = Line1a - Line2b;209 factor = parallel.ScalarProduct(a)/a.Norm();210 if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {211 res = Line2b;212 Log() << Verbose(1) << "Lines conincide." << endl;213 return res;214 }215 }216 Log() << Verbose(1) << "Lines are parallel." << endl;217 res.Zero();218 throw LinearDependenceException(__FILE__,__LINE__);219 }220 221 // obtain s222 double s;223 Vector temp1, temp2;224 temp1 = c;225 temp1.VectorProduct(b);226 temp2 = a;227 temp2.VectorProduct(b);228 Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;229 if (fabs(temp2.NormSquared()) > MYEPSILON)230 s = temp1.ScalarProduct(temp2)/temp2.NormSquared();231 else232 s = 0.;233 Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;234 235 // construct intersection236 res = a;237 res.Scale(s);238 res += Line1a;239 Log() << Verbose(1) << "Intersection is at " << res << "." << endl;240 241 return res;242 };
Note:
See TracChangeset
for help on using the changeset viewer.
