diff --git a/Test/CMakeLists.txt b/Test/CMakeLists.txt deleted file mode 100644 index 7c24083..0000000 --- a/Test/CMakeLists.txt +++ /dev/null @@ -1,99 +0,0 @@ -# @file CMakeLists.txt -# @Author Benjamin Bercovici (bebe0705@colorado.edu) -# @date January, 2017 -# @brief CMake listing enabling compilation of YORP computation example - - - -################################################################################ -# -# User-defined paths -# Should be checked for consistency -# Before running 'cmake ..' in build dir -# -################################################################################ - -# OMP-friendly GCC Major version (will only be used on Mac OS systems) -set (OMP_FRIENDLY_GCC_MAJOR_VERSION 7) - -# Path to OMP-Friendly GCC Compiler (will only be used on Mac OS systems) -set (OMP_FRIENDLY_GCC_PATH /usr/local/Cellar/gcc/7.2.0/bin/) - -# Path to SbgatCoreConfig.cmake (not used if installed at a standard location) -set(YORPLIB_LOC "/home/bebe0705/libs/local/lib/cmake/YORPLib") - -################################################################################ -# -# -# The following should normally not require any modification -# Unless new files are added to the build tree -# -# -################################################################################ - - -if (EXISTS /home/bebe0705/.am_fortuna) - set(IS_FORTUNA ON) - set(RBK_LOC "/home/bebe0705/libs/local/lib/cmake/RigidBodyKinematics") - # set(SBGAT_LOC "/home/bebe0705/libs/local/lib/cmake/SbgatCore") - set(ASPEN_LOC "/home/bebe0705/libs/local/lib/cmake/ASPEN") - -else() - set(IS_FORTUNA OFF) - message("-- This is not Fortuna") - -endif() - - -# Building procedure -get_filename_component(dirName ${CMAKE_CURRENT_SOURCE_DIR} NAME) -set(EXE_NAME ${dirName} CACHE STRING "Name of executable to be created.") - - -if (${IS_FORTUNA}) - set(CMAKE_C_COMPILER "/usr/local/bin/gcc" CACHE STRING "C Compiler" FORCE) - set(CMAKE_CXX_COMPILER "/usr/local/bin/g++" CACHE STRING "C++ Compiler" FORCE) -else() - if(APPLE) - if(EXISTS ${OMP_FRIENDLY_GCC_PATH}) # Else running on a MAC. Will attempt to switch compiler to get the OMP-friendly GCC 6.3.0_1 from Homebrew - message("Switching to OMP-friendly GCC ") - set(CMAKE_C_COMPILER ${OMP_FRIENDLY_GCC_PATH}gcc-${OMP_FRIENDLY_GCC_MAJOR_VERSION} CACHE STRING "C Compiler" FORCE) - set(CMAKE_CXX_COMPILER ${OMP_FRIENDLY_GCC_PATH}g++-${OMP_FRIENDLY_GCC_MAJOR_VERSION} CACHE STRING "C++ Compiler" FORCE) - - else()# no OMP-compliant compiler was found on this mac. - message(FATAL_ERROR "No OMP-compliant compiler was found on this Mac.") - endif() - else() # Running on Linux. Will switch back to compiler in /usr/local/bin - message("Switching to /usr/local/gcc ") - set(CMAKE_C_COMPILER "/usr/local/bin/gcc" CACHE STRING "C Compiler" FORCE) - set(CMAKE_CXX_COMPILER "/usr/local/bin/g++" CACHE STRING "C++ Compiler" FORCE) - endif() -endif() - -project(${EXE_NAME}) - -# Specify the version used -if (${CMAKE_MAJOR_VERSION} LESS 3) - message(FATAL_ERROR " You are running an outdated version of CMake") -endif() - -cmake_minimum_required(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}.0) -set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/source/cmake) - -add_definitions(-Wall -O2 ) -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++17") - -# Find OpenMP -find_package(OpenMP REQUIRED) -set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") - -# Find YORPLib -find_package(YORPLib REQUIRED PATHS ${YORPLIB_LOC}) -include_directories(${YORPLIB_INCLUDE_HEADER}) - -# Add source files in root directory -add_executable(${EXE_NAME} main.cpp) - -target_link_libraries(${EXE_NAME} ${YORPLIB_LIBRARY}) - diff --git a/Test/input_use1.txt b/Test/input_use1.txt deleted file mode 100644 index c790089..0000000 --- a/Test/input_use1.txt +++ /dev/null @@ -1,9 +0,0 @@ -../../../ASPEN_GUI_LESS/resources/shape_models/cube.obj -40 -0.05 0.0 -1.0 -1.0 -2 -3 -5 -./ diff --git a/Test/main.cpp b/Test/main.cpp deleted file mode 100644 index c836e8d..0000000 --- a/Test/main.cpp +++ /dev/null @@ -1,139 +0,0 @@ -// -// main.cpp -// Spacecraft_SRP -// -// Created by Jay McMahon on 5/19/14. -// Copyright (c) 2014 Jay McMahon. All rights reserved. -// - -#include - -#include - -#include -using namespace std; - -int main(int argc, const char * argv[]) -{ - - auto start = std::chrono::system_clock::now(); - - // Inputs have two choices, either a file was input on the command line, or there will be prompts on the screen - // to start, just do the file option - // When implment screen option, write the inputs to a file for re-use - string objFileName, opticalFileName, outputFileBaseName, inFileName; - int numVox, opticalFlag; - double rho, spec; - double lambdaDel; - double deltaDel; - double MaxFourier; - int howManyBounces; - int numrefine; - - if (argc == 1) { - // Screen option used if no extra input arguments - cout << "Enter full path to input file:\n"; - getline (cin, inFileName); - } else if (argc==2) { - // File option used if 1 extra input argument (the file name) - // make this into a sub-function in the future - - // Open the file - inFileName = argv[1]; - - } else { - // Error! - cerr << "This program takes zero or one inputs only\n"; - exit(1); - } - - ifstream inFile(inFileName); - string linestr; - - // First line is .obj file name - getline(inFile, objFileName); - - // Second line is number of voxels per axis - getline(inFile, linestr); - numVox = stoi(linestr); - - // Third line contains either rho and s, or a file name to read from - getline(inFile, linestr); - - // split line based on white space - stringstream ss(linestr); - string item; - vector elems; - while (getline(ss, item, ' ')) { - elems.push_back(item); - } - if (elems.size()==2) { - // Two elements = rho then s - rho = stod(elems[0]); - spec = stod(elems[1]); - opticalFlag = 0; - } else { - // One element = optical property file name - opticalFileName = elems[0]; - opticalFlag = 1; - } - - // Fourth line contains lambdaDel - getline(inFile, linestr); - - lambdaDel = stod(linestr); - - // Fifth line contains deltaDel - getline(inFile, linestr); - deltaDel = stod(linestr); - - // Sixth line contains MaxFourier - getline(inFile, linestr); - MaxFourier = stoi(linestr); - - // Seventh line contains howManyBounces - getline(inFile, linestr); - howManyBounces = stoi(linestr); - - // Eigth line contains numrefine - getline(inFile, linestr); - numrefine = stoi(linestr); - - // Ninth line contains output file base name - getline(inFile, outputFileBaseName); - - - inFile.close(); - - // Construct the Body object with the appropriate method based on opticalFlag - Body targetObj; - if (opticalFlag == 0) { - Body inputObj(objFileName,rho,spec); - targetObj = inputObj; - } else if (opticalFlag == 1) { - Body inputObj(objFileName,opticalFileName); - targetObj = inputObj; - - } - - - // Determine limits for voxel grid - double xmax = 1.01 * targetObj.getMaxDim(1); - double ymax = 1.01 * targetObj.getMaxDim(2); - double zmax = 1.01 * targetObj.getMaxDim(3); - - // Fill out Voxel Grid - targetObj.setVoxelGrid(xmax, ymax, zmax, numVox); - - // Compute and print SRP Fourier coefficients - SRPModel targetSRP(lambdaDel, deltaDel, MaxFourier, &targetObj, howManyBounces, numrefine); - targetSRP.writeSRPCoeffsFile(outputFileBaseName, 2*(90.0/deltaDel) + 1); - - // stop timer and report - auto end = std::chrono::system_clock::now(); - std::chrono::duration elapsed_secs = end - start; - cout << "Time elapsed = " << elapsed_secs.count() << " s" << "\n"; - - return 0; -} - diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt deleted file mode 100644 index bc1eb59..0000000 --- a/lib/CMakeLists.txt +++ /dev/null @@ -1,139 +0,0 @@ -# MIT License - -# Copyright (c) 2014 Jay McMahon - -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: - -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. - -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. -# - -# @file CMakeLists.txt -# @Author Benjamin Bercovici (bebe0705@colorado.edu) -# @date March, 2018 -# @brief CMake listing enabling compilation and installation of the YORPLib library - -################################################################################ -# -# User-defined paths -# Should be checked for consistency -# Before running 'cmake ..' in build dir -# -################################################################################ - -# OMP-friendly GCC Major version (will only be used on Mac OS systems) -set (OMP_FRIENDLY_GCC_MAJOR_VERSION 7) - -# Path to OMP-Friendly GCC Compiler (will only be used on Mac OS systems) -set (OMP_FRIENDLY_GCC_PATH /usr/local/Cellar/gcc/7.2.0/bin/) - -################################################################################ -# -# -# The following should normally not require any modification -# Unless new files are added to the build tree -# -# -################################################################################ - -if (EXISTS /home/bebe0705/.am_fortuna) - set(IS_FORTUNA ON) - message("-- This is Fortuna") - -else() - set(IS_FORTUNA OFF) - message("-- This is not Fortuna") -endif() - - -# Building procedure -get_filename_component(dirName ${CMAKE_CURRENT_SOURCE_DIR} NAME) -set(LIB_NAME YORPLib) - - -if (${IS_FORTUNA}) - set(CMAKE_C_COMPILER "/usr/local/bin/gcc" CACHE STRING "C Compiler" FORCE) - set(CMAKE_CXX_COMPILER "/usr/local/bin/g++" CACHE STRING "C++ Compiler" FORCE) -else() - if(APPLE) - if(EXISTS ${OMP_FRIENDLY_GCC_PATH}) # Else running on a MAC. Will attempt to switch compiler to get the OMP-friendly GCC 6.3.0_1 from Homebrew - message("Switching to OMP-friendly GCC ") - set(CMAKE_C_COMPILER ${OMP_FRIENDLY_GCC_PATH}gcc-${OMP_FRIENDLY_GCC_MAJOR_VERSION} CACHE STRING "C Compiler" FORCE) - set(CMAKE_CXX_COMPILER ${OMP_FRIENDLY_GCC_PATH}g++-${OMP_FRIENDLY_GCC_MAJOR_VERSION} CACHE STRING "C++ Compiler" FORCE) - - else()# no OMP-compliant compiler was found on this mac. - message(FATAL_ERROR "No OMP-compliant compiler was found on this Mac.") - endif() - else() # Running on Linux. Will switch back to compiler in /usr/local/bin - message("Switching to /usr/local/gcc ") - set(CMAKE_C_COMPILER "/usr/local/bin/gcc" CACHE STRING "C Compiler" FORCE) - set(CMAKE_CXX_COMPILER "/usr/local/bin/g++" CACHE STRING "C++ Compiler" FORCE) - endif() -endif() - -project(${LIB_NAME}) - -# Specify the version used -if (${CMAKE_MAJOR_VERSION} LESS 3) - message(FATAL_ERROR " You are running an outdated version of CMake") -endif() - -cmake_minimum_required(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}.0) -set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/source/cmake) - -# Compiler flags -add_definitions(-Wall -O2 ) - -# Enable C++17 -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++17") - -# Include directories -include_directories(include) - -# Find OpenMP -find_package(OpenMP REQUIRED) -set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") - -# Add source files in root directory -add_library(${LIB_NAME} -SHARED -source/Body.cpp -source/Facet.cpp -source/FCoeffs.cpp -source/Shadow.cpp -source/SRPModel.cpp -source/VoxelGrid.cpp -) - -# Installing -if (NOT ${IS_FORTUNA}) - install (TARGETS ${LIB_NAME} DESTINATION /usr/local/lib/) - install (FILES ${PROJECT_SOURCE_DIR}/YORPLibConfig.cmake DESTINATION /usr/local/lib/cmake/YORPLib/) - install (DIRECTORY ${PROJECT_SOURCE_DIR}/include/ DESTINATION /usr/local/include/YORPLib FILES_MATCHING PATTERN "*.h") -else() - install (TARGETS ${LIB_NAME} DESTINATION /home/bebe0705/libs/local/lib) - install (FILES ${PROJECT_SOURCE_DIR}/YORPLibConfig.cmake DESTINATION /home/bebe0705/libs/local/lib/cmake/YORPLib/) - install (DIRECTORY ${PROJECT_SOURCE_DIR}/include/ DESTINATION /home/bebe0705/libs/local/include/YORPLib FILES_MATCHING PATTERN "*.h") -endif() - - - - - - - - diff --git a/lib/YORPLibConfig.cmake b/lib/YORPLibConfig.cmake deleted file mode 100644 index 28c65ca..0000000 --- a/lib/YORPLibConfig.cmake +++ /dev/null @@ -1,51 +0,0 @@ -# MIT License - -# Copyright (c) 2014 Jay McMahon - -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: - -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. - -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. -# - - -if (EXISTS /home/bebe0705/.am_fortuna) - set(IS_FORTUNA ON) - message("-- This is Fortuna") - -else() - set(IS_FORTUNA OFF) -endif() - - -if(${IS_FORTUNA}) - set(YORPLIB_INCLUDE_HEADER /home/bebe0705/libs/local/include/YORPLib/) - set(YORPLIB_LIBRARY /home/bebe0705/libs/local/lib/libYORPLib.so) - -else() - set(YORPLIB_INCLUDE_HEADER /usr/local/include/YORPLib/) - - if (APPLE) - set(YORPLIB_LIBRARY /usr/local/lib/libYORPLib.dylib) - elseif(UNIX AND NOT APPLE) - set(YORPLIB_LIBRARY /usr/local/lib/libYORPLib.so) - else() - message(FATAL_ERROR "Unsupported platform") - endif() -endif() - - -message("-- Found YORPLib: " ${YORPLIB_LIBRARY}) diff --git a/lib/build/.gitignore b/lib/build/.gitignore deleted file mode 100644 index aff905f..0000000 --- a/lib/build/.gitignore +++ /dev/null @@ -1,2 +0,0 @@ -* -*/ diff --git a/lib/include/Body.h b/lib/include/Body.h deleted file mode 100644 index 8ab45c3..0000000 --- a/lib/include/Body.h +++ /dev/null @@ -1,100 +0,0 @@ - -/** MIT License - -Copyright (c) 2014 Jay McMahon - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. -*/ - - - -// -// Body.h -// Spacecraft_SRP -// -// Created by Jay McMahon on 5/19/14. -// Copyright (c) 2014 Jay McMahon. All rights reserved. -// - -#ifndef __Spacecraft_SRP__Body__ -#define __Spacecraft_SRP__Body__ - -//#include -//#include -#include -//#include "Facet.h" -#include "VoxelGrid.h" - -//using namespace std; - -class Body -{ -public: - Body(); - Body(string bodyFile); - Body(string bodyFile, double rho, double spec); - - // Added by Benjamin Bercovici - Body(std::vector > vertices, - std::vector > facets, - double rho, double spec); - - Body(string bodyFile, string optFile); - int getNumFacets(); - int getNumVerts(); -// void getFacetInfo(int fnum, Facet* fac); - double* getFacetNormal(int fnum); - double* getFacetRc(int fnum); - double getFacetArea(int fnum); - int* getNeighbors(int fnum); - vector* getInView(int fnum); - vector* getInViewRc(int fnum); - vector* getF2V(int fnum); - int getVinVox(int voxnum); - int getFinVox(int voxnum); - Facet getFacet( int fnum); - void setNeighbors(); - void setView(); - void setViewRc(); - void setViewGrid(); - void setFacets(double rho, double s); - void setFacetsVec(vector* rho, vector* s); - void setF2V(); - void setVoxelGrid(double xmax_in, double ymax_in, double zmax_in, int N); - int ray_intersection(double* pt, double* uhat, double* hit_pnt, vector* extraFacets); - double getMaxDimVoxGrid(); -// void getOpticalFourier(int fnum, double& rho, double& s, double& a2, double* nn, double* Ar1, double* Ar2, double* Ar3, double* Ar1_2, double* Ar2_2, double* Ar3_2); - void getOpticalFourier(int fnum, double& rho, double& s, double& a2, double nn[3][3], double Ar1[3], double Ar2[3][3], double Ar3[3][3], double Ar1_2[3], double Ar2_2[3][3], double Ar3_2[3][3]); - double getMaxDim(int axis); -private: - int numFacets, numVerts; - bool viewCheck; // true if setView has been run, false otherwise - vector neighbors; - vector < vector* > inView; - vector < vector* > inViewRc; - vector < vector* > inViewGrid; - vector facetList; // contains the facet numbers as referenced IN THE FILE, so the first vertex = 1... - vector vertices; - vector facetData; - vector < vector > f2vlist; // list of facets associate with each vertex (index 1 is associated with vertex - VoxelGrid bodyVox; - double dot ( double* a, double* b); -}; - -#endif /* defined(__Spacecraft_SRP__Body__) */ diff --git a/lib/include/FCoeffs.h b/lib/include/FCoeffs.h deleted file mode 100644 index 5185d6b..0000000 --- a/lib/include/FCoeffs.h +++ /dev/null @@ -1,100 +0,0 @@ -/** MIT License - -Copyright (c) 2014 Jay McMahon - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. -*/ - - - -// -// FCoeffs.h -// Spacecraft_SRP -// -// Created by Jay McMahon on 5/27/14. -// Copyright (c) 2014 Jay McMahon. All rights reserved. -// - -#ifndef __Spacecraft_SRP__FourierCoeffs__ -#define __Spacecraft_SRP__FourierCoeffs__ - -//#include -//#include -//#include -//#include "Facet.h" -//#include "VoxelGrid.h" -//#include "Body.h" -#include "Shadow.h" - -class FCoeffs -{ -public: - FCoeffs(); - FCoeffs(int order); - void computeCoeffs(const double delta_s, Body* Target, const vector < vector >* low, const vector < vector >* high, const double lamDel, const int numbounces, const vector < vector >* fmult, const vector < vector >* tmult, const vector < vector < double > >* latSpec, const vector < vector < double > >* longSpec, const vector < vector < vector < double > > >* fSpec, const vector < vector < vector < double > > >* tSpec); - void computeYear(); - void setOrder(int order); - double* getA(); - double* getB(); - double* getC(); - double* getD(); - -private: - int order; - double A[3]; - double B[3]; - double Ayear[3]; - double Byear[3]; - double C[3]; - double D[3]; - double Cyear[3]; - double Dyear[3]; - double compute_uc1n(double lbnd, double ubnd, int nint); - double compute_uc2n(double lbnd, double ubnd, int nint); - double compute_uc3n(double lbnd, double ubnd, int nint); - double compute_uuc11n(double lbnd, double ubnd, int nint); - double compute_uuc12n(double lbnd, double ubnd, int nint); - double compute_uuc22n(double lbnd, double ubnd, int nint); - double compute_us1n(double lbnd, double ubnd, int nint); - double compute_us2n(double lbnd, double ubnd, int nint); - double compute_us3n(double lbnd, double ubnd, int nint); - double compute_uus11n(double lbnd, double ubnd, int nint); - double compute_uus12n(double lbnd, double ubnd, int nint); - double compute_uus22n(double lbnd, double ubnd, int nint); - double compute_uc1n_shad(double a, double b, int nint, double Ha, double dH); - double compute_uc2n_shad(double a, double b, int nint, double Ha, double dH); - double compute_uc3n_shad(double a, double b, int nint, double Ha, double dH); - double compute_uuc11n_shad(double a, double b, int nint, double Ha, double dH); - double compute_uuc12n_shad(double a, double b, int nint, double Ha, double dH); - double compute_uuc22n_shad(double a, double b, int nint, double Ha, double dH); - double compute_us1n_shad(double a, double b, int nint, double Ha, double dH); - double compute_us2n_shad(double a, double b, int nint, double Ha, double dH); - double compute_us3n_shad(double a, double b, int nint, double Ha, double dH); - double compute_uus11n_shad(double a, double b, int nint, double Ha, double dH); - double compute_uus12n_shad(double a, double b, int nint, double Ha, double dH); - double compute_uus22n_shad(double a, double b, int nint, double Ha, double dH); - double compute_u1_shad(double a, double b, double Ha, double dH); - double compute_u2_shad(double a, double b, double Ha, double dH); - double compute_u3_shad(double a, double b, double Ha, double dH); - double compute_uu11_shad(double a, double b, double Ha, double dH); - double compute_uu12_shad(double a, double b, double Ha, double dH); - double compute_uu22_shad(double a, double b, double Ha, double dH); -}; - -#endif /* defined(__Spacecraft_SRP__FourierCoeffs__) */ diff --git a/lib/include/Facet.h b/lib/include/Facet.h deleted file mode 100644 index ad00e19..0000000 --- a/lib/include/Facet.h +++ /dev/null @@ -1,80 +0,0 @@ -/** MIT License - -Copyright (c) 2014 Jay McMahon - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. -*/ - - - -// -// Facet.h -// Spacecraft_SRP -// -// Created by Jay McMahon on 5/19/14. -// Copyright (c) 2014 Jay McMahon. All rights reserved. -// - -#ifndef __Spacecraft_SRP__Facet__ -#define __Spacecraft_SRP__Facet__ - -#include -#include - -using namespace std; - -class Facet -{ -public: - Facet(); - Facet(int* verts, vector* vertices, double rho, double s); - void updateRhoS(double rhoSet, double sSet); - void computeOptical(); - int getVert(int vnum); - double getRho(); - double getS(); - double getArea(); - double* getRc(); - double* getNormal(); - double geta2(); - void getOpticalFourier(double nnOut[3][3], double Ar1Out[3], double Ar2Out[3][3], double Ar3Out[3][3], double Ar1_2Out[3], double Ar2_2Out[3][3], double Ar3_2Out[3][3]); -// void getOpticalFourier(double* nnOut, double* Ar1Out, double* Ar2Out, double* Ar3Out, double* Ar1_2Out, double* Ar2_2Out, double* Ar3_2Out); -private: - int vert_nums[3]; // index number of first element of each vertex (vertex 1 spans indices 0-2 in the vector) - //double v1[3], v2[3], v3[3]; - double r_c[3]; - double normal[3]; - double rho; - double s; - double area; -// Below are the properties needed for the Lambertian BRDF; consider packaging these into some sort of BRDF class so that it will be easier in the future to switch out BRDFs -// double B; // Assumed = 2/3 for now; just doing Lambertian to start - double a2; - double nn[3][3]; // = normals(:,j)*normals(:,j)'; - double Ar1[3]; // = -areas(j)./pi.*normals(:,j); - double Ar2[3][3]; // = -areas(j)./pi.*nn; - double Ar3[3][3]; // = -areas(j)./pi.*(2.*nn - eye(3)); - double Ar1_2[3]; // = Ar1./2; - double Ar2_2[3][3]; // = Ar2./2; - double Ar3_2[3][3]; // = Ar3./2; - // Private member functions - void computeGeometry(double* v1, double* v2, double* v3, double* normal, double* r_c, double* area); -}; - -#endif /* defined(__Spacecraft_SRP__Facet__) */ diff --git a/lib/include/SRPModel.h b/lib/include/SRPModel.h deleted file mode 100644 index 6c25129..0000000 --- a/lib/include/SRPModel.h +++ /dev/null @@ -1,94 +0,0 @@ -/** MIT License - -Copyright (c) 2014 Jay McMahon - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. -*/ - - -// -// SRPModel.h -// Spacecraft_SRP -// -// Created by Jay McMahon on 5/27/14. -// Copyright (c) 2014 Jay McMahon. All rights reserved. -// - -#ifndef __Spacecraft_SRP__SRPModel__ -#define __Spacecraft_SRP__SRPModel__ - -//#include -//#include -//#include -#include -#include "FCoeffs.h" -//#include "Facet.h" -//#include "VoxelGrid.h" -//#include "Body.h" -//#include "Shadow.h" - -//using namespace std; - -class SRPModel -{ -public: - SRPModel(); - SRPModel(double lambdaDel, double deltaDel, int MaxFourier, Body* bodyExist, int bounces, int numRefine); - void writeSRPCoeffs(int deltaWrite); - void writeSRPCoeffsFile(string outputFileBaseName, int numDelta); -private: - Body Spacecraft; - vector deltaSunList; - vector lambdaSunList; - int numBounces, numGrid, numShadRefine; // if these are 0, then not using this logic (do this to start) - int FourierOrder; - vector shadowing; // first index corresponds to lambda_s, each shadow has the value for all facets - vector < vector > shadowBoundsLow; // 1st facet, 2nd bound case - vector < vector > shadowBoundsHigh; - vector < vector > riseLambda; // outer index corresponds to delta_s, inner to each facet - vector < vector > setLambda; // outer index corresponds to delta_s, inner to each facet - vector < vector > Coefficients; // outer index corresponds to delta_s, inner to Fourier order - void computeMulti(); - vector < vector < double > > fMult; - vector < vector < double > > tMult; - vector < vector < vector < double > > > fSpec; // facet, specular case num, force vector - vector < vector < vector < double > > > tSpec; // facet, specular case num, torque vector - vector < vector < double > > latSpec; // facet, solar latitude for each specular case num - vector < vector < double > > longSpec; // facet, solar longitude for each specular case num - void computeRiseSet(const double* normal, const double latitude, double* rise, double* set); - double dot ( double* a, double* b); - void cross(double* a, double* b, double* c); - void specularMulti(const int dsnum, const int lsnum); - void compute_fi(const double* uhat, Facet hitF, double* fi); - void computeShadowBounds(const int dsindex, const int numf, const int lsnum, const double lambdaDel); -// double I0_1[3]; // = zeros(3,1); -// double I0_2[3][3]; // = zeros(3); -// vector Ic_1; // = cell(fourier_max,1); -// vector Is_1; // = cell(fourier_max,1); -// vector Ic_2; // = cell(fourier_max,1); -// vector Is_2; // = cell(fourier_max,1); -/* for k = 1:fourier_max - Ic_1{k} = zeros(3,1); - Is_1{k} = zeros(3,1); - Ic_2{k} = zeros(3); - Is_2{k} = zeros(3); - end */ -}; - -#endif /* defined(__Spacecraft_SRP__SRPModel__) */ diff --git a/lib/include/Shadow.h b/lib/include/Shadow.h deleted file mode 100644 index 495654f..0000000 --- a/lib/include/Shadow.h +++ /dev/null @@ -1,57 +0,0 @@ -/** MIT License - -Copyright (c) 2014 Jay McMahon - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. -*/ - - -// -// Shadow.h -// Spacecraft_SRP -// -// Created by Jay McMahon on 5/27/14. -// Copyright (c) 2014 Jay McMahon. All rights reserved. -// - -#ifndef __Spacecraft_SRP__Shadow__ -#define __Spacecraft_SRP__Shadow__ - -//#include -//#include -#include "Body.h" - -//using namespace std; - -class Shadow -{ -public: - Shadow(); - Shadow(int numfacets, double delta, double lambda); - void ComputeShadowing(Body* Target, vector* riseLam, vector* setLam); - void ComputeShadowingGrid(Body* Target, vector* riseLam, vector* setLam); // add grid defining inputs - int getPcntLit(int fnum); -private: - double deltaSun; - double lambdaSun; - vector percentUnShadowed; // Need percent of area that sees the Sun for coeff calculations - double dot ( double* a, double* b); -}; - -#endif /* defined(__Spacecraft_SRP__Shadow__) */ diff --git a/lib/include/VoxelGrid.h b/lib/include/VoxelGrid.h deleted file mode 100644 index 9297c59..0000000 --- a/lib/include/VoxelGrid.h +++ /dev/null @@ -1,82 +0,0 @@ -/** MIT License - -Copyright (c) 2014 Jay McMahon - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. -*/ - - -// -// VoxelGrid.h -// Spacecraft_SRP -// -// Created by Jay McMahon on 5/19/14. -// Copyright (c) 2014 Jay McMahon. All rights reserved. -// - -#ifndef __Spacecraft_SRP__VoxelGrid__ -#define __Spacecraft_SRP__VoxelGrid__ - -//#include -//#include -#include "Facet.h" - -//using namespace std; - -class VoxelGrid -{ -public: - VoxelGrid(); - VoxelGrid(double xmax_in, double ymax_in, double zmax_in, int N); - void fillGrid(vector* facets, vector* vertices, int numVerts, vector< vector >* f2vlist, vector* fdata); - int ray_intersect_voxel(double* pt, double* uhat, vector* facets, vector* vertices, vector * facetData, double* hit_pnt, vector* extraFacets); - void find_pt_in_voxel(double* pt, vector* voxel_list); - int test_intersection(double* pt, double* uhat, vector* facets, vector* vertices, vector * facetData, double*hit_pnt, int current_voxel, vector* extraFacets); - int getNumFinVox(int VoxNum); - int getNumVinVox(int VoxNum); - double getMaxDim(); - void getVoxInd(int voxnum, int* x, int* y, int* z); - bool intersectTriVox(int voxNum, double* v1, double* v2, double* v3, double* nhat); -private: - vector xmax; - vector xmin; - vector ymax; - vector ymin; - vector zmax; - vector zmin; - vector posx; - vector negx; - vector posy; - vector negy; - vector posz; - vector negz; - vector fnum; - vector vnum; - vector < vector > flist; - vector < vector > vlist; - double dx, dy, dz; - int N; - void inverse ( double* answer, double* x, double* y, double* n); - bool in_triangle_3D( double *A, double *B, double *C, double *P, double *nhat); - double dot ( double* a, double* b); - bool satAABBTri(double rBox, double* axis, double* v1, double* v2, double* v3); - void cross(double* a, double* b, double* c); -}; - -#endif /* defined(__Spacecraft_SRP__VoxelGrid__) */ diff --git a/lib/source/Body.cpp b/lib/source/Body.cpp deleted file mode 100644 index 075c577..0000000 --- a/lib/source/Body.cpp +++ /dev/null @@ -1,737 +0,0 @@ -/** MIT License - -Copyright (c) 2014 Jay McMahon - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. -*/ - - -// -// Body.cpp -// Spacecraft_SRP -// -// Created by Jay McMahon on 5/19/14. -// Copyright (c) 2014 Jay McMahon. All rights reserved. -// - -#include "Body.h" -#include -#include -#include -#include -#include -#include -#include - -using namespace std; - -Body::Body(){ - // cerr << "Need inputs to define Body\n"; // change this to throw an error - - return; -} - -Body::Body(string bodyFile, double rho, double spec) -{ - ifstream inFile ( bodyFile ); - numFacets = 0; - numVerts = 0; - - if ( !inFile.is_open() ) { - // The file could not be opened - cerr << "Requested body file doesn't exist\n"; // change this to throw an error - } - else { - - - // Safely use the file stream - string linestr, item; - vector elems; - - // Assume that the input file is a .obj file, so the format has lines that start with "v" or "f" - // get a line - while (getline(inFile, linestr)) { - // split line based on white space - stringstream ss(linestr); - while (getline(ss, item, ' ')) { - elems.push_back(item); - } - - if (elems[0]=="v") { - vertices.push_back(stod(elems[1])); - vertices.push_back(stod(elems[2])); - vertices.push_back(stod(elems[3])); - numVerts++; - } else if (elems[0]=="f") { - facetList.push_back(stoi(elems[1])); - facetList.push_back(stoi(elems[2])); - facetList.push_back(stoi(elems[3])); - numFacets++; - } - elems.clear(); - } - } - inFile.close(); - - // Call other member functions to fill out the rest of the object - facetData.resize(numFacets); - - setFacets(rho, spec); - neighbors.resize(3*numFacets); - - setNeighbors(); - inView.resize(numFacets); - - setView(); - f2vlist.resize(numVerts); - setF2V(); - -} - - -Body::Body(std::vector > vertices, - std::vector > facets, double rho, double spec){ - - numFacets = facets.size(); - numVerts = vertices.size(); - - - for (int i = 0; i < vertices.size(); ++i){ - this -> vertices.push_back(vertices[i][0]); - this -> vertices.push_back(vertices[i][1]); - this -> vertices.push_back(vertices[i][2]); - } - - for (int i = 0; i < facets.size(); ++i){ - this -> facetList.push_back(facets[i][0] + 1); - this -> facetList.push_back(facets[i][1] + 1); - this -> facetList.push_back(facets[i][2] + 1); - } - - - // Call other member functions to fill out the rest of the object - facetData.resize(numFacets); - - setFacets(rho, spec); - neighbors.resize(3*numFacets); - - setNeighbors(); - inView.resize(numFacets); - - setView(); - f2vlist.resize(numVerts); - setF2V(); - -} - - - - - - -Body::Body(string bodyFile) -{ - ifstream inFile ( bodyFile ); - numFacets = 0; - numVerts = 0; - - if ( !inFile.is_open() ) { - // The file could not be opened - cerr << "Requested body file doesn't exist\n"; // change this to throw an error - } - else { - // Safely use the file stream - string linestr, item; - vector elems; - - // Assume that the input file is a .obj file, so the format has lines that start with "v" or "f" - // get a line - while (getline(inFile, linestr)) { - // split line based on white space - stringstream ss(linestr); - while (getline(ss, item, ' ')) { - elems.push_back(item); - } - if (elems[0]=="v") { - vertices.push_back(stod(elems[1])); - vertices.push_back(stod(elems[2])); - vertices.push_back(stod(elems[3])); - numVerts++; - } else if (elems[0]=="f") { - facetList.push_back(stoi(elems[1])); - facetList.push_back(stoi(elems[2])); - facetList.push_back(stoi(elems[3])); - std::cout << facetList[numFacets - 3] << " " << facetList[numFacets - 2] << " " << facetList[numFacets - 1] << std::endl; - numFacets++; - } - elems.clear(); - } - } - inFile.close(); - - // Call other member functions to fill out the rest of the object - facetData.resize(numFacets); - setFacets(0.0, 0.0); - neighbors.resize(3*numFacets); - setNeighbors(); - inView.resize(numFacets); - setView(); - setF2V(); -} - -Body::Body(string bodyFile, string optFile) -{ - ifstream inFile ( bodyFile ); - ifstream inFile2 ( optFile ); - numFacets = 0; - numVerts = 0; - - if ( !inFile.is_open() ) { - // The file could not be opened - cerr << "Requested body file doesn't exist\n"; // change this to throw an error - } - else { - // Safely use the file stream - string linestr, item; - vector elems; - - // Assume that the input file is a .obj file, so the format has lines that start with "v" or "f" - // get a line - while (getline(inFile, linestr)) { - // split line based on white space - stringstream ss(linestr); - while (getline(ss, item, ' ')) { - elems.push_back(item); - } - if (elems[0]=="v") { - vertices.push_back(stod(elems[1])); - vertices.push_back(stod(elems[2])); - vertices.push_back(stod(elems[3])); - numVerts++; - } else if (elems[0]=="f") { - facetList.push_back(stoi(elems[1])); - facetList.push_back(stoi(elems[2])); - facetList.push_back(stoi(elems[3])); - numFacets++; - } - elems.clear(); - } - } - inFile.close(); - - vector rhoVec(numFacets); - vector specVec(numFacets); - - if ( !inFile2.is_open() ) { - // The file could not be opened - cerr << "Requested optical parameter file doesn't exist\n"; // change this to throw an error - } - else { - // Safely use the file stream - string linestr, item; - vector elems; - int counter = 0; - - // This file has one line for every facet, each line has two double values: reflectivity (rho) and the percent of that which is specular (s) - while (getline(inFile2, linestr)) { - // split line based on white space - stringstream ss(linestr); - while (getline(ss, item, ' ')) { - elems.push_back(item); - } - rhoVec[counter] = stod(elems[0]); - specVec[counter] = stod(elems[1]); - counter++; - elems.clear(); - } - } - inFile2.close(); - - // Call other member functions to fill out the rest of the object - facetData.resize(numFacets); - setFacetsVec(&rhoVec, &specVec); - neighbors.resize(3*numFacets); - setNeighbors(); - inView.resize(numFacets); - setView(); - f2vlist.resize(numVerts); - setF2V(); -} - -void Body::setNeighbors() -{ - int i, j; - int x, y, z, a, b, c; - int xc, yc, zc; - vector numneigh(numFacets,0); - - for (i=0;i* rho, vector* s) -{ - int ii; - int verts[3]; - - for (ii=0; ii* rowP = new vector; - //row.resize((numFacets-1)/2); - //rowP->resize(2000); // 4530 per is max? - inView[ii] = rowP; - } - - // Check for facets that COULD see each other - // Test relative to all vertices to account for all alignments; break as soon as a true statement is found -// int vCount; - for (ii=0; ii* rowP = new vector; - rowP->resize((numFacets-ii+5)/5); - inView[ii] = rowP; */ -// vCount = 0; - -// for (jj=ii+1; jj 1.0e-10) && (dot(&rtemp[0], facetData[jj].getNormal()) < -1.0e-10)) { // JWM 2/10/16 took out >=/<= to get rid of planar facets seeing eachother - // JWM 2/24/16 normalized rtemp; set limit so vertex has to be > 1e-10 radians above horizon to count as in view -// try { - //if (vCount >= inView[ii]->size()) { - // inView[ii]->resize(inView[ii]->size() + 2000); - //cout << "size inVew[" << ii << "] = " << inView[ii]->size() << "\n"; - //} - //inView[ii]->at(vCount) = jj; - inView[ii]->push_back(jj); -// vCount++; - //inView[ii]->push_back(jj); - //inView[jj]->push_back(ii); -// } catch (...) { -// cout << "size inVew[" << ii << "] = " << inView[ii]->size() << "\n"; -// cout << "size inVew[" << jj << "] = " << inView[jj]->size() << "\n"; -// cout << "jay\n"; -// } - flag = true; - break; - } - } - if (flag==true) { - break; - } - } - -// if (vCount > 0) { -// break; -// } - } - - - - //inView[ii]->shrink_to_fit(); - - /*if (ii==77){ - for (mm=0; mmsize(); mm++) { - cout << inView[ii]->at(mm) << "\n"; - } - cout << "jay\n"; - }*/ - - } - - -// for (ii=0; iishrink_to_fit(); -// } - - viewCheck = true; - -} - -void Body::setViewRc() -{ - // Check for intervening facets based only on center-to-center view - int ii, jj, kk, test; - double hit_pnt[3]; - double uhat[3], vmag, pt[3]; - double* rc1; - double* rc2; - vector dummy; - - // If for some reason the view hasn't been set yet, do so now - if (!viewCheck) { - setView(); - } - - // Initialize inViewRc as a copy of inView - inViewRc.resize(numFacets); - for (ii=0; ii* rowP = new vector; - inViewRc[ii] = rowP; - for (jj=0; jj<(*inView[ii]).size(); jj++){ - inViewRc[ii]->push_back((*inView[ii])[jj]); - } - } - - // Check every center-to-center ray trace to see if it intersects the facet that is "inView" - for (ii=0; iiempty()) { - rc1 = facetData[ii].getRc(); - for (jj=0; jjsize(); jj++) { - // skip if this connection already checked because these are symmetric - if (ii > inViewRc[ii]->at(jj)) { - continue; - } - - // set the test direction from facet ii to the facet in inViewRc[ii][jj] - rc2 = facetData[inViewRc[ii]->at(jj)].getRc(); - for (kk=0; kk<3; kk++) { - uhat[kk] = rc2[kk] - rc1[kk]; - } - vmag = sqrt(dot(&uhat[0], &uhat[0])); // JWM 2/10/16 moved normalization up here so the horizon checking would actually work! - for (kk=0; kk<3; kk++) { - uhat[kk] = uhat[kk]/vmag; - } - - // ray tracing doesn't detect if it goes below either facets' horizon - use previous test for this - // JWM 2/24/16 modified to have to be > 1e-10 radians above horizon to count as in view - if ((dot(&uhat[0], facetData[ii].getNormal()) < -1.0e-10) || (dot(&uhat[0], facetData[inViewRc[ii]->at(jj)].getNormal()) > 1.0e-10)) { - // This means the center connecting line goes through the facets, so erase - // Erase the mirror inViewRc entry also - inViewRc[inViewRc[ii]->at(jj)]->erase(find(inViewRc[inViewRc[ii]->at(jj)]->begin(), inViewRc[inViewRc[ii]->at(jj)]->end(), ii)); - inViewRc[ii]->erase(inViewRc[ii]->begin()+jj); - jj--; - } else { - for (kk=0; kk<3; kk++) { - pt[kk] = rc1[kk] + 0.001*vmag*uhat[kk]; // move starting point 0.1% of the distance along the connecting line so the ray tracing check doesn't just return the starting facet - } - test = bodyVox.ray_intersect_voxel(&pt[0], &uhat[0], &facetList, &vertices, &facetData, &hit_pnt[0], &dummy); - if (inViewRc[ii]->at(jj) != test) { - // This means we intersected a different facet before the targeted one, so need to erase - // Erase the mirror inViewRc entry also - // Note that we don't check for multiple facets because if it hits the facet it is supposed to, it won't hit multiple facets - inViewRc[inViewRc[ii]->at(jj)]->erase(find(inViewRc[inViewRc[ii]->at(jj)]->begin(), inViewRc[inViewRc[ii]->at(jj)]->end(), ii)); - inViewRc[ii]->erase(inViewRc[ii]->begin()+jj); - jj--; - } - } - } - } - } -} - -void Body::setViewGrid() -{ - // Using a grid requires extra outputs beyond just inViewGrid, which says what facets can see each other at all. - // Also need to output how much they can see of one another, and the center of pressure offsets from the original center of facet. - - // Check for intervening facets based on vertex-to-vertex views, and compute visible portions - int ii, jj, kk, test, mm, pp, vCount; - double hit_pnt[3]; - double uhat[3], vmag, pt[3]; - double* rc1; - double* rc2; - vector dummy; - - // If for some reason the view hasn't been set yet, do so now - if (!viewCheck) { - setView(); - } - - // If grid hasn't been initialized do so now - - inViewGrid = inView; - - // Check every vertex-to-vertex ray trace to see if it intersects the facet that is "inView" - for (ii=0; iiempty()) { - - for (jj=0; jjsize(); jj++) { - // skip if this connection already checked because these are symmetric - if (ii > inViewGrid[ii]->at(jj)) { - continue; - } - vCount = 0; - for (mm=0; mm<3; mm++) { - rc1 = &vertices[facetData[ii].getVert(mm)]; - for (pp=0; pp<3; pp++) { - // set the test direction from facet ii to the facet in inViewGrid[ii][jj] - rc2 = &vertices[facetData[inViewGrid[ii]->at(jj)].getVert(pp)]; - for (kk=0; kk<3; kk++) { - uhat[kk] = rc2[kk] - rc1[kk]; - } - vmag = sqrt(dot(&uhat[0], &uhat[0])); - for (kk=0; kk<3; kk++) { - uhat[kk] = uhat[kk]/vmag; - pt[kk] = rc1[kk] + 0.001*vmag*uhat[kk]; // move starting point 0.1% of the distance along the connecting line so the ray tracing check doesn't just return the starting facet - } - test = bodyVox.ray_intersect_voxel(&pt[0], &uhat[0], &facetList, &vertices, &facetData, &hit_pnt[0], &dummy); - if (test == inViewGrid[ii]->at(jj)) { - ++vCount; - } - } - } - - - if (vCount == 0) { - // This means we intersected a different facet before the targeted one, so need to erase - // Erase the mirror inViewGrid entry also - inViewGrid[inViewGrid[ii]->at(jj)]->erase(find(inViewGrid[inViewGrid[ii]->at(jj)]->begin(), inViewGrid[inViewGrid[ii]->at(jj)]->end(), ii)); - inViewGrid[ii]->erase(inViewGrid[ii]->begin()+jj); - } else if (vCount < 9) { - // This means the view is partially blocked - compute grid outputs - - } - } - } - } -} - - -void Body::setF2V() -{ - vector::iterator iter_f; - vector < vector >::iterator iter_f2v; - vector row; - int ii = 1; // starts at 1 because facetlist has vertex numbers starting at 1 - int jj; - - for (iter_f2v = f2vlist.begin(); iter_f2v != f2vlist.end(); ++iter_f2v) { - *iter_f2v = row; - iter_f = facetList.begin(); - while (iter_f != facetList.end()) { - iter_f = find(iter_f, facetList.end(), ii); - if (iter_f != facetList.end()) { - jj = iter_f - facetList.begin(); - iter_f2v->push_back(jj/3); // will give the vertex number, starting at 0 - ++iter_f; - } - } - ii++; - } - -} - - -void Body::setVoxelGrid(double xmax_in, double ymax_in, double zmax_in, int N) -{ - VoxelGrid bodyVoxTemp(xmax_in, ymax_in, zmax_in, N); - bodyVox = bodyVoxTemp; - bodyVox.fillGrid(&facetList, &vertices, numVerts, &f2vlist, &facetData); -} - -int Body::ray_intersection(double* pt, double* uhat, double* hit_pnt, vector* extraFacets) -{ - int test_hit; - test_hit = bodyVox.ray_intersect_voxel(pt, uhat, &facetList, &vertices, &facetData, hit_pnt, extraFacets); - return test_hit; -} - -double Body::dot ( double* a, double* b) -{ - return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; - -} - -int Body::getNumFacets() -{ - return numFacets; -} - -int Body::getNumVerts() -{ - return numVerts; -} - -/*void Body::getFacetInfo(int fnum, Facet* fac) -{ - cout << &facetData.at(fnum); - fac = &facetData.at(fnum); - return; -}*/ - -Facet Body::getFacet( int fnum) -{ - return facetData[fnum]; -} - -double* Body::getFacetNormal(int fnum) -{ - return facetData[fnum].getNormal(); -} - -double* Body::getFacetRc(int fnum) -{ - return facetData[fnum].getRc(); -} - -double Body::getFacetArea(int fnum) -{ - return facetData[fnum].getArea(); -} - -int Body::getVinVox(int voxnum) -{ - return bodyVox.getNumVinVox(voxnum); -} - -int Body::getFinVox(int voxnum) -{ - return bodyVox.getNumFinVox(voxnum); -} - -int* Body::getNeighbors(int fnum) -{ - return &neighbors[3*fnum]; -} - -vector* Body::getInView(int fnum) -{ - return inView[fnum]; -} - -vector* Body::getInViewRc(int fnum) -{ - return inViewRc[fnum]; -} - -vector* Body::getF2V(int fnum) -{ - return &f2vlist[fnum]; -} - -double Body::getMaxDimVoxGrid() -{ - return bodyVox.getMaxDim(); -} - -void Body::getOpticalFourier(int fnum, double& rho, double& s, double& a2, double nn[3][3], double Ar1[3], double Ar2[3][3], double Ar3[3][3], double Ar1_2[3], double Ar2_2[3][3], double Ar3_2[3][3]) -{ -// double temp; - a2 = facetData[fnum].geta2(); -// a2 = &temp; -// temp = facetData[fnum].getRho(); -// rho = &temp; - rho = facetData[fnum].getRho(); - s = facetData[fnum].getS(); -// s = &temp; - - facetData[fnum].getOpticalFourier(nn, Ar1, Ar2, Ar3, Ar1_2, Ar2_2, Ar3_2); - - return; - -} - -double Body::getMaxDim(int axis) -{ - // axis == 1 is x axis - // axis == 2 is y axis - // axis == 3 is z axis - - int offset; - - if (axis == 1) { - offset = 0; - } else if (axis == 2) { - offset = 1; - } else if (axis == 3) { - offset = 2; - } - - vector data; - for (int ii = offset; ii<3*numVerts; ii+=3) { - data.push_back(abs(vertices[ii])); - } - - return *max_element(data.begin(), data.end()); - -} \ No newline at end of file diff --git a/lib/source/FCoeffs.cpp b/lib/source/FCoeffs.cpp deleted file mode 100644 index 920060d..0000000 --- a/lib/source/FCoeffs.cpp +++ /dev/null @@ -1,794 +0,0 @@ -/** MIT License - -Copyright (c) 2014 Jay McMahon - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. -*/ - - -// -// FCoeffs.cpp -// Spacecraft_SRP -// -// Created by Jay McMahon on 5/27/14. -// Copyright (c) 2014 Jay McMahon. All rights reserved. -// - -#include "FCoeffs.h" -#include "Body.h" -#include -#include - -FCoeffs::FCoeffs() -{ - FCoeffs(1); -} - -FCoeffs::FCoeffs(int orderIn) -{ - order = orderIn; - for (int ii=0; ii<3; ii++) { - A[ii] = 0.0; - B[ii] = 0.0; - Ayear[ii] = 0.0; - Byear[ii] = 0.0; - C[ii] = 0.0; - D[ii] = 0.0; - Cyear[ii] = 0.0; - Dyear[ii] = 0.0; - } -} - -void FCoeffs::computeCoeffs(const double delta_s, Body* Target, const vector < vector >* low, const vector < vector >* high, const double lamDel, const int numbounces, const vector < vector >* fmult, const vector < vector >* tmult, const vector < vector < double > >* latSpec, const vector < vector < double > >* longSpec, const vector < vector < vector < double > > >* fSpec, const vector < vector < vector < double > > >* tSpec) -{ - int numf = Target->getNumFacets(); - double cdels = cos(delta_s); - double sdels = sin(delta_s); - double cdels2 = cdels*cdels; - double sdels2 = sdels*sdels; - double csdels = cdels*sdels; - double u1, u2, u3, uu11, uu12, uu22; - double I0_1[3], I0_2[3][3]; - int jj, kk, ll, bound_count; - double rho, s, a2; - double nn[3][3], Ar1[3], Ar2[3][3], Ar3[3][3]; - double Ar1_2[3], Ar2_2[3][3], Ar3_2[3][3]; - double nhat[3], rc[3]; - double* nhat_p; - double* rc_p; -// double* rho_p; - double A0[3], A1[3], A2[3], temp[3]; - double B0[3], B1[3], B2[3], temp2[3]; - double uc1n, uc2n, uc3n, us1n, us2n, us3n; - double uuc11n, uuc12n, uuc22n, uus11n, uus12n, uus22n; - double Ic_1[3], Is_1[3], Ic_2[3][3], Is_2[3][3]; - double A01m, A1m, B1m; - vector ::const_iterator iter_d; - int specIndex; - - for (jj=0; jj<3; jj++) { - A[jj] = 0.0; - B[jj] = 0.0; - C[jj] = 0.0; - D[jj] = 0.0; - } - - for (int ii=0; iigetOpticalFourier(ii, rho, s, a2, nn, Ar1, Ar2, Ar3, Ar1_2, Ar2_2, Ar3_2); - nhat_p = Target->getFacetNormal(ii); - nhat[0] = *nhat_p; - nhat[1] = *(nhat_p + 1); - nhat[2] = *(nhat_p + 2); - rc_p = Target->getFacetRc(ii); - rc[0] = *rc_p; - rc[1] = *(rc_p+1); - rc[2] = *(rc_p+2); - // get number of segments - bound_count = (*low)[ii].size(); - double dHall = 1.0/lamDel; - if (order==0) { - // Compute integrals analytically over each unshadowed segment - u1 = 0.0; - u2 = 0.0; - u3 = 0.0; - uu11 = 0.0; - uu12 = 0.0; - uu22 = 0.0; - for (kk = 0; kk 0.0) { - u1 += compute_u1_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], 0.0, dHall); - u2 += compute_u2_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], 0.0, dHall); - u3 += compute_u3_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], 0.0, dHall); - uu11 += compute_uu11_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], 0.0, dHall); - uu12 += compute_uu12_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], 0.0, dHall); - uu22 += compute_uu22_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], 0.0, dHall); - } - // Compute the partial shadowing after the high boundary - if ((*high)[ii][kk] < 2.0*M_PI) { - u1 += compute_u1_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, 1.0, -dHall); - u2 += compute_u2_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, 1.0, -dHall); - u3 += compute_u3_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, 1.0, -dHall); - uu11 += compute_uu11_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, 1.0, -dHall); - uu12 += compute_uu12_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, 1.0, -dHall); - uu22 += compute_uu22_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, 1.0, -dHall); - } - } - - // Compute zero order integrals - I0_1[0] = cdels*u1; - I0_1[1] = cdels*u2; - I0_1[2] = sdels*u3; - - I0_2[0][0] = cdels2*uu11; - I0_2[1][0] = cdels2*uu12; - I0_2[2][0] = csdels*u1; - I0_2[0][1] = I0_2[1][0]; - I0_2[1][1] = cdels2*uu22; - I0_2[2][1] = csdels*u2; - I0_2[0][2] = I0_2[2][0]; - I0_2[1][2] = I0_2[2][1]; - I0_2[2][2] = sdels2*u3; - - // Compute Total Coefficients - A01m = 0.0; - for (jj=0; jj<3; jj++) { - A0[jj] = 0.0; - A1[jj] = 0.0; - A2[jj] = 0.0; - for (kk=0; kk<3; kk++) { - A0[jj] += I0_2[jj][kk]*Ar1_2[kk]; - A1[jj] += Ar2_2[jj][kk]*I0_1[kk]; - for (ll=0; ll<3; ll++) { - A2[jj] += Ar3_2[jj][kk]*I0_2[kk][ll]*nhat[ll]; - } - if ((numbounces>0) && (jj==0)) { - A01m += Ar1_2[kk]*I0_1[kk]; - } - } - temp[jj] = A0[jj] + a2*A1[jj] + rho*s*A2[jj]; - A[jj] += temp[jj]; - if (numbounces>0) { - A[jj] += A01m*(*fmult)[ii][jj]; - } - } - C[0] += rc[1]*temp[2] - rc[2]*temp[1]; - C[1] += rc[2]*temp[0] - rc[0]*temp[2]; - C[2] += rc[0]*temp[1] - rc[1]*temp[0]; - if (numbounces>0) { - for (jj=0; jj<3; jj++) { - C[jj] += A01m*(*tmult)[ii][jj]; - } - } - } else { - // Compute integrals analytically over each unshadowed segment - if (order == 1) { - u1 = 0.0; - u2 = 0.0; - uu11 = 0.0; - uu12 = 0.0; - uu22 = 0.0; - for (kk = 0; kk 0.0) { - u1 += compute_u1_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], 0.0, dHall); - u2 += compute_u2_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], 0.0, dHall); - uu11 += compute_uu11_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], 0.0, dHall); - uu12 += compute_uu12_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], 0.0, dHall); - uu22 += compute_uu22_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], 0.0, dHall); - } - // Compute the partial shadowing after the high boundary - if ((*high)[ii][kk] < 2.0*M_PI) { - u1 += compute_u1_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, 1.0, -dHall); - u2 += compute_u2_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, 1.0, -dHall); - uu11 += compute_uu11_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, 1.0, -dHall); - uu12 += compute_uu12_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, 1.0, -dHall); - uu22 += compute_uu22_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, 1.0, -dHall); - } - } - uc1n = uu11; - uc2n = uu12; - uc3n = u1; - us1n = uu12; - us2n = uu22; - us3n = u2; - } else { - uc1n = 0.0; - uc2n = 0.0; - uc3n = 0.0; - us1n = 0.0; - us2n = 0.0; - us3n = 0.0; - for (kk = 0; kk 0.0) { - uc1n += compute_uc1n_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], order, 0.0, dHall); - uc2n += compute_uc2n_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], order, 0.0, dHall); - uc3n += compute_uc3n_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], order, 0.0, dHall); - us1n += compute_us1n_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], order, 0.0, dHall); - us2n += compute_us2n_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], order, 0.0, dHall); - us3n += compute_us3n_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], order, 0.0, dHall); - } - // Compute the partial shadowing after the high boundary - if ((*high)[ii][kk] < 2.0*M_PI) { - uc1n += compute_uc1n_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, order, 1.0, -dHall); - uc2n += compute_uc2n_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, order, 1.0, -dHall); - uc3n += compute_uc3n_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, order, 1.0, -dHall); - us1n += compute_us1n_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, order, 1.0, -dHall); - us2n += compute_us2n_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, order, 1.0, -dHall); - us3n += compute_us3n_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, order, 1.0, -dHall); - } - } - } - uuc11n = 0.0; - uuc12n = 0.0; - uuc22n = 0.0; - uus11n = 0.0; - uus12n = 0.0; - uus22n = 0.0; - for (kk = 0; kk 0.0) { - uuc11n += compute_uuc11n_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], order, 0.0, dHall); - uuc12n += compute_uuc12n_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], order, 0.0, dHall); - uuc22n += compute_uuc22n_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], order, 0.0, dHall); - uus11n += compute_uus11n_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], order, 0.0, dHall); - uus12n += compute_uus12n_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], order, 0.0, dHall); - uus22n += compute_uus22n_shad((*low)[ii][kk] - lamDel, (*low)[ii][kk], order, 0.0, dHall); - } - // Compute the partial shadowing after the high boundary - if ((*high)[ii][kk] < 2.0*M_PI) { - uuc11n += compute_uuc11n_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, order, 1.0, -dHall); - uuc12n += compute_uuc12n_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, order, 1.0, -dHall); - uuc22n += compute_uuc22n_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, order, 1.0, -dHall); - uus11n += compute_uus11n_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, order, 1.0, -dHall); - uus12n += compute_uus12n_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, order, 1.0, -dHall); - uus22n += compute_uus22n_shad((*high)[ii][kk], (*high)[ii][kk]+lamDel, order, 1.0, -dHall); - } - } - - // Compute nth order integrals - // Cosine 1 Integral - Ic_1[0] = cdels*uc1n; - Ic_1[1] = cdels*uc2n; - Ic_1[2] = sdels*uc3n; - - // Sine 1 Integral - Is_1[0] = cdels*us1n; - Is_1[1] = cdels*us2n; - Is_1[2] = sdels*us3n; - - // Cosine 2 Integral - Ic_2[0][0] = cdels2*uuc11n; - Ic_2[1][0] = cdels2*uuc12n; - Ic_2[2][0] = csdels*uc1n; - Ic_2[0][1] = Ic_2[1][0]; - Ic_2[1][1] = cdels2*uuc22n; - Ic_2[2][1] = csdels*uc2n; - Ic_2[0][2] = Ic_2[2][0]; - Ic_2[1][2] = Ic_2[2][1]; - Ic_2[2][2] = sdels2*uc3n; - - // Sine 2 Integral - Is_2[0][0] = cdels2*uus11n; - Is_2[1][0] = cdels2*uus12n; - Is_2[2][0] = csdels*us1n; - Is_2[0][1] = Is_2[1][0]; - Is_2[1][1] = cdels2*uus22n; - Is_2[2][1] = csdels*us2n; - Is_2[0][2] = Is_2[2][0]; - Is_2[1][2] = Is_2[2][1]; - Is_2[2][2] = sdels2*us3n; - - // Compute Total Coefficients - A1m = 0.0; - B1m = 0.0; - for (jj=0; jj<3; jj++) { - A0[jj] = 0.0; - A1[jj] = 0.0; - A2[jj] = 0.0; - B0[jj] = 0.0; - B1[jj] = 0.0; - B2[jj] = 0.0; - for (kk=0; kk<3; kk++) { - A0[jj] += Ic_2[jj][kk]*Ar1[kk]; - A1[jj] += Ar2[jj][kk]*Ic_1[kk]; - B0[jj] += Is_2[jj][kk]*Ar1[kk]; - B1[jj] += Ar2[jj][kk]*Is_1[kk]; - for (ll=0; ll<3; ll++) { - A2[jj] += Ar3[jj][kk]*Ic_2[kk][ll]*nhat[ll]; - B2[jj] += Ar3[jj][kk]*Is_2[kk][ll]*nhat[ll]; - } - if ((numbounces>0) && (jj==0)) { - A1m += Ar1[kk]*Ic_1[kk]; - B1m += Ar1[kk]*Is_1[kk]; - } - } - temp[jj] = A0[jj] + a2*A1[jj] + rho*s*A2[jj]; - A[jj] += temp[jj]; - - temp2[jj] = B0[jj] + a2*B1[jj] + rho*s*B2[jj]; - B[jj] += temp2[jj]; - - if (numbounces>0) { - A[jj] += A1m*(*fmult)[ii][jj]; - B[jj] += B1m*(*fmult)[ii][jj]; - } - } - C[0] += rc[1]*temp[2] - rc[2]*temp[1]; - C[1] += rc[2]*temp[0] - rc[0]*temp[2]; - C[2] += rc[0]*temp[1] - rc[1]*temp[0]; - - D[0] += rc[1]*temp2[2] - rc[2]*temp2[1]; - D[1] += rc[2]*temp2[0] - rc[0]*temp2[2]; - D[2] += rc[0]*temp2[1] - rc[1]*temp2[0]; - - if (numbounces>0) { - for (jj=0; jj<3; jj++) { - C[jj] += A1m*(*tmult)[ii][jj]; - D[jj] += B1m*(*tmult)[ii][jj]; - } - } - - } - } - - // add specular reflections - if ((numbounces > 0) && (!(*latSpec)[ii].empty())) { - - iter_d = find((*latSpec)[ii].begin(),(*latSpec)[ii].end(), delta_s); - - while (iter_d != (*latSpec)[ii].end()) { - // get associated index - specIndex = iter_d - (*latSpec)[ii].begin(); - - // Add contribution from matching latitutde specular reflections to current order coefficients - if (order == 0) { - for (jj=0; jj<3; jj++) { - A[jj] += (*fSpec)[ii][specIndex][jj]/(2*M_PI); - C[jj] += (*tSpec)[ii][specIndex][jj]/(2*M_PI); - } - } else { - for (jj=0; jj<3; jj++) { - A[jj] += (*fSpec)[ii][specIndex][jj] * cos(order * (*longSpec)[ii][specIndex])/M_PI; - B[jj] += (*fSpec)[ii][specIndex][jj] * sin(order * (*longSpec)[ii][specIndex])/M_PI; - C[jj] += (*tSpec)[ii][specIndex][jj] * cos(order * (*longSpec)[ii][specIndex])/M_PI; - D[jj] += (*tSpec)[ii][specIndex][jj] * sin(order * (*longSpec)[ii][specIndex])/M_PI; - } - } - - // update iterator for next search - iter_d = find(iter_d+1,(*latSpec)[ii].end(), delta_s); - } - } - } -} - -double FCoeffs::compute_uc1n(double lbnd, double ubnd, int nint) -{ - double n = (double)nint; - double fu = (1.0/(n*n-1.0))*(-cos(n*ubnd)*sin(ubnd) + n*cos(ubnd)*sin(n*ubnd)); - double fl = (1.0/(n*n-1.0))*(-cos(n*lbnd)*sin(lbnd) + n*cos(lbnd)*sin(n*lbnd)); - - return fu - fl; - -} - -double FCoeffs::compute_uc2n(double lbnd, double ubnd, int nint) -{ - double n = (double)nint; - double fu = (1.0/(n*n-1.0))*(cos(n*ubnd)*cos(ubnd) + n*sin(ubnd)*sin(n*ubnd)); - double fl = (1.0/(n*n-1.0))*(cos(n*lbnd)*cos(lbnd) + n*sin(lbnd)*sin(n*lbnd)); - - return fu - fl; - -} - -double FCoeffs::compute_uc3n(double lbnd, double ubnd, int nint) -{ - double n = (double)nint; - double fu = (1.0/n)*(sin(n*ubnd)); - double fl = (1.0/n)*(sin(n*lbnd)); - - return fu - fl; - -} - -double FCoeffs::compute_uuc11n(double lbnd, double ubnd, int nint) -{ - double n = (double)nint; - double fu, fl; - if (nint != 2) { - fu = (1.0/4.0)*((sin((n-2.0)*ubnd)/(n-2.0)) + (2.0*sin(n*ubnd)/n) + (sin((n+2.0)*ubnd)/(n+2.0))); - fl = (1.0/4.0)*((sin((n-2.0)*lbnd)/(n-2.0)) + (2.0*sin(n*lbnd)/n) + (sin((n+2.0)*lbnd)/(n+2.0))); - } else { - fu = (ubnd/4.0) + (sin(2.0*ubnd)/4.0) + (sin(4.0*ubnd)/16.0); - fl = (lbnd/4.0) + (sin(2.0*lbnd)/4.0) + (sin(4.0*lbnd)/16.0); - } - - return fu - fl; - -} - -double FCoeffs::compute_uuc12n(double lbnd, double ubnd, int nint) -{ - double n = (double)nint; - double fu, fl; - if (nint != 2) { - fu = (1.0/4.0)*((cos((n-2.0)*ubnd)/(n-2.0)) - (cos((n+2.0)*ubnd)/(n+2.0))); - fl = (1.0/4.0)*((cos((n-2.0)*lbnd)/(n-2.0)) - (cos((n+2.0)*lbnd)/(n+2.0))); - } else { - fu = -cos(4.0*ubnd)/16.0; - fl = -cos(4.0*lbnd)/16.0; - } - - return fu - fl; - -} - -double FCoeffs::compute_uuc22n(double lbnd, double ubnd, int nint) -{ - double n = (double)nint; - double fu, fl; - - if (nint != 2) { - fu = (1.0/4.0)*((-sin((n-2.0)*ubnd)/(n-2.0)) + (2.0*sin(n*ubnd)/n) - (sin((n+2.0)*ubnd)/(n+2.0))); - fl = (1.0/4.0)*((-sin((n-2.0)*lbnd)/(n-2.0)) + (2.0*sin(n*lbnd)/n) - (sin((n+2.0)*lbnd)/(n+2.0))); - } else { - fu = (-ubnd/4.0) + (sin(2.0*ubnd)/4.0) - (sin(4.0*ubnd)/16.0); - fl = (-lbnd/4.0) + (sin(2.0*lbnd)/4.0) - (sin(4.0*lbnd)/16.0); - } - - return fu - fl; - -} - -double FCoeffs::compute_us1n(double lbnd, double ubnd, int nint) -{ - double n = (double)nint; - double fu = (1.0/(1.0-n*n))*(sin(n*ubnd)*sin(ubnd) + n*cos(ubnd)*cos(n*ubnd)); - double fl = (1.0/(1.0-n*n))*(sin(n*lbnd)*sin(lbnd) + n*cos(lbnd)*cos(n*lbnd)); - - return fu - fl; - -} - -double FCoeffs::compute_us2n(double lbnd, double ubnd, int nint) -{ - double n = (double)nint; - double fu = (1.0/(n*n-1.0))*(-n*cos(n*ubnd)*sin(ubnd) + cos(ubnd)*sin(n*ubnd)); - double fl = (1.0/(n*n-1.0))*(-n*cos(n*lbnd)*sin(lbnd) + cos(lbnd)*sin(n*lbnd)); - - return fu - fl; - -} - -double FCoeffs::compute_us3n(double lbnd, double ubnd, int nint) -{ - double n = (double)nint; - double fu = (1.0/n)*(-cos(n*ubnd)); - double fl = (1.0/n)*(-cos(n*lbnd)); - - return fu - fl; - -} - -double FCoeffs::compute_uus11n(double lbnd, double ubnd, int nint) -{ - double n = (double)nint; - double fu, fl; - - if (nint != 2) { - fu = (1.0/4.0)*((-cos((n-2.0)*ubnd)/(n-2.0)) - (2.0*cos(n*ubnd)/n) - (cos((n+2.0)*ubnd)/(n+2.0))); - fl = (1.0/4.0)*((-cos((n-2.0)*lbnd)/(n-2.0)) - (2.0*cos(n*lbnd)/n) - (cos((n+2.0)*lbnd)/(n+2.0))); - } else { - fu = -cos(ubnd)*cos(ubnd)*cos(ubnd)*cos(ubnd)/2.0; - fl = -cos(lbnd)*cos(lbnd)*cos(lbnd)*cos(lbnd)/2.0; - } - - return fu - fl; - -} - -double FCoeffs::compute_uus12n(double lbnd, double ubnd, int nint) -{ - double n = (double)nint; - double fu, fl; - - if (n != 2) { - fu = (1.0/4.0)*((sin((n-2.0)*ubnd)/(n-2.0)) - (sin((n+2.0)*ubnd)/(n+2.0))); - fl = (1.0/4.0)*((sin((n-2.0)*lbnd)/(n-2.0)) - (sin((n+2.0)*lbnd)/(n+2.0))); - } else { - fu = (ubnd/4.0) - (sin(4.0*ubnd)/16.0); - fl = (lbnd/4.0) - (sin(4.0*lbnd)/16.0); - } - - return fu - fl; - -} - -double FCoeffs::compute_uus22n(double lbnd, double ubnd, int nint) -{ - double n = (double)nint; - double fu, fl; - - if (n != 2) { - fu = (1.0/4.0)*((cos((n-2.0)*ubnd)/(n-2.0)) - (2.0*cos(n*ubnd)/n) + (cos((n+2.0)*ubnd)/(n+2.0))); - fl = (1.0/4.0)*((cos((n-2.0)*lbnd)/(n-2.0)) - (2.0*cos(n*lbnd)/n) + (cos((n+2.0)*lbnd)/(n+2.0))); - } else { - fu = sin(ubnd)*sin(ubnd)*sin(ubnd)*sin(ubnd)/2.0; - fl = sin(lbnd)*sin(lbnd)*sin(lbnd)*sin(lbnd)/2.0; - } - - return fu - fl; - -} - -double FCoeffs::compute_uc1n_shad(double a, double b, int nint, double Ha, double dH) -{ - double n = (double)nint; - double ans; - - if (nint == 1) { - ans = ((-(dH*cos(2*a)) + dH*cos(2*b) + 2*((a - b)*(a*dH - b*dH - 2*Ha) - Ha*sin(2*a) + (-(a*dH) + b*dH + Ha)*sin(2*b))))/8.0; - } else { - ans = ((-(dH*cos(a*(-1 + n))) + dH*cos(b*(-1 + n)) + ((-1 + n)*(-(dH*(-1 + n)*cos(a*(1 + n))) + 2*Ha*(1 + n)*(cos(a*n)*sin(a) - n*cos(a)*sin(a*n))))/pow(1 + n,2) + ((-1 + n)*(-(sin(b)*(2*(-(a*dH) + b*dH + Ha)*(1 + n)*cos(b*n) + dH*(-1 + n)*sin(b*n))) + cos(b)*(dH*(-1 + n)*cos(b*n) + 2*(-(a*dH) + b*dH + Ha)*n*(1 + n)*sin(b*n))))/pow(1 + n,2)))/(2.*pow(-1 + n,2)); - } - - return ans; -} - -double FCoeffs::compute_uc2n_shad(double a, double b, int nint, double Ha, double dH) -{ - double n = (double)nint; - double ans; - - if (nint == 1) { - ans = ((2*Ha*cos(2*a) + 2*(a*dH - b*dH - Ha)*cos(2*b) + dH*(-sin(2*a) + sin(2*b))))/8.0; - } else { - ans = ((-((Ha*cos(a*(-1 + n)))/(-1 + n)) + ((-(a*dH) + b*dH + Ha)*cos(b*(-1 + n)))/(-1 + n) + (Ha*cos(a*(1 + n)))/(1 + n) - ((-(a*dH) + b*dH + Ha)*cos(b*(1 + n)))/(1 + n) + (dH*sin(a*(-1 + n)))/pow(-1 + n,2) - (dH*sin(b*(-1 + n)))/pow(-1 + n,2) - (dH*sin(a*(1 + n)))/pow(1 + n,2) + (dH*sin(b*(1 + n)))/pow(1 + n,2)))/2.0; - } - - return ans; -} - -double FCoeffs::compute_uc3n_shad(double a, double b, int nint, double Ha, double dH) -{ - double n = (double)nint; - double ans; - - ans = ((-(dH*cos(a*n)) + dH*cos(b*n) - Ha*n*sin(a*n) + (-(a*dH) + b*dH + Ha)*n*sin(b*n)))/pow(n,2); - - return ans; -} - -double FCoeffs::compute_us1n_shad(double a, double b, int nint, double Ha, double dH) -{ - double n = (double)nint; - double ans; - - if (nint == 1) { - ans = ((2*Ha*cos(2*a) + 2*(a*dH - b*dH - Ha)*cos(2*b) + dH*(-sin(2*a) + sin(2*b))))/8.0; - } else { - ans = (((Ha*cos(a*(-1 + n)))/(-1 + n) - ((-(a*dH) + b*dH + Ha)*cos(b*(-1 + n)))/(-1 + n) + (Ha*cos(a*(1 + n)))/(1 + n) - ((-(a*dH) + b*dH + Ha)*cos(b*(1 + n)))/(1 + n) - (dH*sin(a*(-1 + n)))/pow(-1 + n,2) + (dH*sin(b*(-1 + n)))/pow(-1 + n,2) - (dH*sin(a*(1 + n)))/pow(1 + n,2) + (dH*sin(b*(1 + n)))/pow(1 + n,2)))/2.0; - } - - return ans; -} - -double FCoeffs::compute_us2n_shad(double a, double b, int nint, double Ha, double dH) -{ - double n = (double)nint; - double ans; - - if (nint == 1) { - ans = ((dH*cos(2*a) - dH*cos(2*b) + 2*(Ha*(sin(2*a) - sin(2*b)) + (a - b)*(a*dH - b*dH - 2*Ha + dH*sin(2*b)))))/8.0; - } else { - ans = ((-(dH*cos(a*(-1 + n))) + dH*cos(b*(-1 + n)) + ((-1 + n)*(dH*(-1 + n)*cos(a*(1 + n)) + 2*Ha*(1 + n)*(n*cos(a*n)*sin(a) - cos(a)*sin(a*n))))/pow(1 + n,2) + ((-1 + n)*(-(sin(b)*(2*(-(a*dH) + b*dH + Ha)*n*(1 + n)*cos(b*n) - dH*(-1 + n)*sin(b*n))) + cos(b)*(-(dH*(-1 + n)*cos(b*n)) + 2*(-(a*dH) + b*dH + Ha)*(1 + n)*sin(b*n))))/pow(1 + n,2)))/(2.*pow(-1 + n,2)); - } - - return ans; -} - -double FCoeffs::compute_us3n_shad(double a, double b, int nint, double Ha, double dH) -{ - double n = (double)nint; - double ans; - - ans = ((Ha*n*cos(a*n) + (a*dH - b*dH - Ha)*n*cos(b*n) + dH*(-sin(a*n) + sin(b*n))))/pow(n,2); - - return ans; -} - -double FCoeffs::compute_uuc11n_shad(double a, double b, int nint, double Ha, double dH) -{ - double n = (double)nint; - double ans; - - if (nint == 2) { - ans = ((8*a*(a*dH - 2*Ha) + 8*b*(-2*a*dH + b*dH + 2*Ha) - 8*dH*cos(2*a) - dH*cos(4*a) + 8*dH*cos(2*b) + dH*cos(4*b) - 4*Ha*(4*sin(2*a) + sin(4*a)) + 4*(-(a*dH) + b*dH + Ha)*(4*sin(2*b) + sin(4*b))))/64.0; - } else { - ans = ((-((dH*cos(a*(-2 + n)))/pow(-2 + n,2)) + (dH*cos(b*(-2 + n)))/pow(-2 + n,2) - (2*dH*cos(a*n))/pow(n,2) + (2*dH*cos(b*n))/pow(n,2) - (dH*cos(a*(2 + n)))/pow(2 + n,2) + (dH*cos(b*(2 + n)))/pow(2 + n,2) - (Ha*sin(a*(-2 + n)))/(-2 + n) - (a*dH*sin(b*(-2 + n)))/(-2 + n) + (b*dH*sin(b*(-2 + n)))/(-2 + n) + (Ha*sin(b*(-2 + n)))/(-2 + n) - (2*Ha*sin(a*n))/n - (2*a*dH*sin(b*n))/n + (2*b*dH*sin(b*n))/n + (2*Ha*sin(b*n))/n - (Ha*sin(a*(2 + n)))/(2 + n) - (a*dH*sin(b*(2 + n)))/(2 + n) + (b*dH*sin(b*(2 + n)))/(2 + n) + (Ha*sin(b*(2 + n)))/(2 + n)))/4.0; - } - - return ans; -} - -double FCoeffs::compute_uuc12n_shad(double a, double b, int nint, double Ha, double dH) -{ - double n = (double)nint; - double ans; - - if (nint == 2) { - ans = ((4*Ha*cos(4*a) - 4*(-(a*dH) + b*dH + Ha)*cos(4*b) - dH*sin(4*a) + dH*sin(4*b)))/64.0; - } else { - ans = ((-((Ha*cos(a*(-2 + n)))/(-2 + n)) + ((-(a*dH) + b*dH + Ha)*cos(b*(-2 + n)))/(-2 + n) + (Ha*cos(a*(2 + n)))/(2 + n) - ((-(a*dH) + b*dH + Ha)*cos(b*(2 + n)))/(2 + n) + (dH*sin(a*(-2 + n)))/pow(-2 + n,2) - (dH*sin(b*(-2 + n)))/pow(-2 + n,2) - (dH*sin(a*(2 + n)))/pow(2 + n,2) + (dH*sin(b*(2 + n)))/pow(2 + n,2)))/4.0; - } - - return ans; -} - -double FCoeffs::compute_uuc22n_shad(double a, double b, int nint, double Ha, double dH) -{ - double n = (double)nint; - double ans; - - if (nint == 2) { - ans = ((-8*a*(a*dH - 2*Ha) - 8*b*(-2*a*dH + b*dH + 2*Ha) - 8*dH*cos(2*a) + dH*cos(4*a) + 8*dH*cos(2*b) - dH*cos(4*b) + 4*Ha*(-4*sin(2*a) + sin(4*a)) - 4*(-(a*dH) + b*dH + Ha)*(-4*sin(2*b) + sin(4*b))))/64.0; - } else { - ans = ((-((a*dH*sin(a*(-2 + n)))/(-2 + n)) + (Ha*sin(a*(-2 + n)))/(-2 + n) + (dH*(cos(a*(-2 + n)) + a*(-2 + n)*sin(a*(-2 + n))))/pow(-2 + n,2) + (a*dH*sin(b*(-2 + n)))/(-2 + n) - (Ha*sin(b*(-2 + n)))/(-2 + n) - (dH*(cos(b*(-2 + n)) + b*(-2 + n)*sin(b*(-2 + n))))/pow(-2 + n,2) + (2*a*dH*sin(a*n))/n - (2*Ha*sin(a*n))/n - (2*dH*(cos(a*n) + a*n*sin(a*n)))/pow(n,2) - (2*a*dH*sin(b*n))/n + (2*Ha*sin(b*n))/n + (2*dH*(cos(b*n) + b*n*sin(b*n)))/pow(n,2) - (a*dH*sin(a*(2 + n)))/(2 + n) + (Ha*sin(a*(2 + n)))/(2 + n) + (dH*(cos(a*(2 + n)) + a*(2 + n)*sin(a*(2 + n))))/pow(2 + n,2) + (a*dH*sin(b*(2 + n)))/(2 + n) - (Ha*sin(b*(2 + n)))/(2 + n) - (dH*(cos(b*(2 + n)) + b*(2 + n)*sin(b*(2 + n))))/pow(2 + n,2)))/4.0; - } - - return ans; -} - -double FCoeffs::compute_uus11n_shad(double a, double b, int nint, double Ha, double dH) -{ - double n = (double)nint; - double ans; - - if (nint == 2) { - ans = ((4*Ha*(4*cos(2*a) + cos(4*a) - 4*cos(2*b) - cos(4*b)) + dH*(4*(a - b)*(4*cos(2*b) + cos(4*b)) - 8*sin(2*a) - sin(4*a) + 8*sin(2*b) + sin(4*b))))/64.0; - } else { - ans = (((Ha*cos(a*(-2 + n)))/(-2 + n) - ((-(a*dH) + b*dH + Ha)*cos(b*(-2 + n)))/(-2 + n) + (2*Ha*cos(a*n))/n - (2*(-(a*dH) + b*dH + Ha)*cos(b*n))/n + (Ha*cos(a*(2 + n)))/(2 + n) + (a*dH*cos(b*(2 + n)))/(2 + n) - (b*dH*cos(b*(2 + n)))/(2 + n) - (Ha*cos(b*(2 + n)))/(2 + n) - (dH*sin(a*(-2 + n)))/pow(-2 + n,2) + (dH*sin(b*(-2 + n)))/pow(-2 + n,2) - (2*dH*sin(a*n))/pow(n,2) + (2*dH*sin(b*n))/pow(n,2) - (dH*sin(a*(2 + n)))/pow(2 + n,2) + (dH*sin(b*(2 + n)))/pow(2 + n,2)))/4.0; - } - - return ans; -} - -double FCoeffs::compute_uus12n_shad(double a, double b, int nint, double Ha, double dH) -{ - double n = (double)nint; - double ans; - - if (nint == 2) { - ans = ((dH*cos(4*a) - dH*cos(4*b) + 4*(Ha*(sin(4*a) - sin(4*b)) + (a - b)*(2*a*dH - 2*b*dH - 4*Ha + dH*sin(4*b)))))/64.0; - } else { - ans = ((-(dH*cos(a*(-2 + n))) + dH*cos(b*(-2 + n)) + ((-2 + n)*(dH*(-2 + n)*cos(a*(2 + n)) + 4*Ha*(2 + n)*(n*cos(a)*cos(a*n)*sin(a) - pow(cos(a),2)*sin(a*n) + pow(sin(a),2)*sin(a*n))))/pow(2 + n,2) - ((-2 + n)*(dH*(-2 + n)*cos(b*(2 + n)) + 4*(-(a*dH) + b*dH + Ha)*(2 + n)*(n*cos(b)*cos(b*n)*sin(b) - pow(cos(b),2)*sin(b*n) + pow(sin(b),2)*sin(b*n))))/pow(2 + n,2)))/(4.0*pow(-2 + n,2)); - } - - return ans; -} - -double FCoeffs::compute_uus22n_shad(double a, double b, int nint, double Ha, double dH) -{ - double n = (double)nint; - double ans; - - if (nint == 2) { - ans = ((4*Ha*(4*cos(2*a) - cos(4*a) - 4*cos(2*b) + cos(4*b)) + dH*(-4*(a - b)*(-4*cos(2*b) + cos(4*b)) - 8*sin(2*a) + sin(4*a) + 8*sin(2*b) - sin(4*b))))/64.0; - } else { - ans = ((-((Ha*cos(a*(-2 + n)))/(-2 + n)) + ((-(a*dH) + b*dH + Ha)*cos(b*(-2 + n)))/(-2 + n) + (2*Ha*cos(a*n))/n - (2*(-(a*dH) + b*dH + Ha)*cos(b*n))/n - (Ha*cos(a*(2 + n)))/(2 + n) - (a*dH*cos(b*(2 + n)))/(2 + n) + (b*dH*cos(b*(2 + n)))/(2 + n) + (Ha*cos(b*(2 + n)))/(2 + n) + (dH*sin(a*(-2 + n)))/pow(-2 + n,2) - (dH*sin(b*(-2 + n)))/pow(-2 + n,2) - (2*dH*sin(a*n))/pow(n,2) + (2*dH*sin(b*n))/pow(n,2) + (dH*sin(a*(2 + n)))/pow(2 + n,2) - (dH*sin(b*(2 + n)))/pow(2 + n,2)))/4.0; - } - - return ans; -} - -double FCoeffs::compute_u1_shad(double a, double b, double Ha, double dH) -{ - double ans = -(dH*cos(a)) + dH*cos(b) - Ha*sin(a) + (-(a*dH) + b*dH + Ha)*sin(b); - return ans; - -} - -double FCoeffs::compute_u2_shad(double a, double b, double Ha, double dH) -{ - double ans = Ha*cos(a) + (a*dH - b*dH - Ha)*cos(b) + dH*(-sin(a) + sin(b)); - return ans; -} - -double FCoeffs::compute_u3_shad(double a, double b, double Ha, double dH) -{ - double ans = (a - b)*(a*dH - b*dH - 2*Ha)/2.0; - return ans; -} - -double FCoeffs::compute_uu11_shad(double a, double b, double Ha, double dH) -{ - double ans = (-(dH*cos(2*a)) + dH*cos(2*b) + 2*((a - b)*(a*dH - b*dH - 2*Ha) - Ha*sin(2*a) + (-(a*dH) + b*dH + Ha)*sin(2*b)))/8.0; - return ans; - -} - -double FCoeffs::compute_uu12_shad(double a, double b, double Ha, double dH) -{ - double ans = (2*Ha*cos(2*a) + 2*(a*dH - b*dH - Ha)*cos(2*b) + dH*(-sin(2*a) + sin(2*b)))/8.0; - return ans; -} - -double FCoeffs::compute_uu22_shad(double a, double b, double Ha, double dH) -{ - double ans = (dH*cos(2*a) - dH*cos(2*b) + 2*(Ha*(sin(2*a) - sin(2*b)) + (a - b)*(a*dH - b*dH - 2*Ha + dH*sin(2*b))))/8.0; - return ans; -} - -// Set the order -void FCoeffs::setOrder(int orderIn) { - order = orderIn; -} - -// Get A -double* FCoeffs::getA() -{ - return &A[0]; -} - -// Get B -double* FCoeffs::getB() -{ - return &B[0]; -} - -// Get C -double* FCoeffs::getC() -{ - return &C[0]; -} - -// Get D -double* FCoeffs::getD() -{ - return &D[0]; -} \ No newline at end of file diff --git a/lib/source/Facet.cpp b/lib/source/Facet.cpp deleted file mode 100644 index 57ac3d2..0000000 --- a/lib/source/Facet.cpp +++ /dev/null @@ -1,229 +0,0 @@ -/** MIT License - -Copyright (c) 2014 Jay McMahon - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. -*/ - - - -// -// Facet.cpp -// Spacecraft_SRP -// -// Created by Jay McMahon on 5/19/14. -// Copyright (c) 2014 Jay McMahon. All rights reserved. -// - -#include "Facet.h" -#include -#include - -// Overridden default constructor -Facet::Facet() -{ - vert_nums[0] = 0; - vert_nums[1] = 0; - vert_nums[2] = 0; - rho = 0.0; - s = 0.0; -} - -// Preferred constructor -Facet::Facet(int* verts, vector* vertices, double rhoSet, double sSet) -{ - // Note that verts numbers start at 1, not 0! may want to change this on input... - vert_nums[0] = 3*(verts[0]-1); - vert_nums[1] = 3*(verts[1]-1); - vert_nums[2] = 3*(verts[2]-1); - - rho = rhoSet; - - s = sSet; - - // Get the vertex coordinates - double* v1; - double* v2; - double* v3; - - // std::cout << "before computeGeometry" << std::endl; - - - // std::cout << verts[0] << " " << verts[1] <<" "<< verts[2] << std::endl; - - // std::cout << vert_nums[0] << " " << vert_nums[1] <<" "<< vert_nums[2] << std::endl; - - - - v1 = &vertices->at(vert_nums[0]); - v2 = &vertices->at(vert_nums[1]); - v3 = &vertices->at(vert_nums[2]); - - // std::cout << "computeGeometry" << std::endl; - // Call member functions to fill remaining values - computeGeometry(&v1[0], &v2[0], &v3[0], &normal[0], &r_c[0], &area); - -} - -// Private function to compute geometry -void Facet::computeGeometry(double* v1, double* v2, double* v3, double* normal, double* r_c, double* area) -{ - int ii; - - // Compute r_c - for (ii=0; ii<3; ii++) { - r_c[ii] = (v1[ii]+v2[ii]+v3[ii])/3; - } - - // Compute edge vectors - double e1[3]; - double e2[3]; - double e3[3]; - for (ii=0; ii<3; ii++) { - e1[ii] = v2[ii] - v1[ii]; - e2[ii] = v3[ii] - v2[ii]; - e3[ii] = v1[ii] - v3[ii]; - } - - - // Compute normal - // This could benefit from a cross product and norm function - normal[0] = e1[1]*e2[2] - e1[2]*e2[1]; - normal[1] = e1[2]*e2[0] - e1[0]*e2[2]; - normal[2] = e1[0]*e2[1] - e1[1]*e2[0]; - double nmag; - nmag = sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]); - for (ii=0; ii<3; ii++) { - normal[ii] = normal[ii]/nmag; - } - - // Compute area - double a, b, c, perim_2; - a = sqrt(e1[0]*e1[0] + e1[1]*e1[1] + e1[2]*e1[2]); - b = sqrt(e2[0]*e2[0] + e2[1]*e2[1] + e2[2]*e2[2]); - c = sqrt(e3[0]*e3[0] + e3[1]*e3[1] + e3[2]*e3[2]); - perim_2 = (a+b+c)/2; - - *area = sqrt(perim_2*(perim_2-a)*(perim_2-b)*(perim_2-c)); -} - -// Get one of the vertex numbers -int Facet::getVert(int vnum) -{ - return vert_nums[vnum]; -} - -// Get reflectivity -double Facet::getRho() -{ - return rho; -} - -// Get specularity -double Facet::getS() -{ - return s; -} - -// Get area -double Facet::getArea() -{ - return area; -} - -// Get facet center -double* Facet::getRc() -{ - return &r_c[0]; -} - -// Get normal -double* Facet::getNormal() -{ - return &normal[0]; -} - -// Get a2 parameter -double Facet::geta2() -{ - return a2; -} - -// Set Lambertian model optical properties needed to compute coefficients -void Facet::computeOptical() -{ - // Lambertian - double B = 2.0/3.0; - - a2 = B*rho*(1.0-s) + B*(1.0-rho); - - int ii, jj; - - for (ii=0; ii<3; ii++) { - Ar1[ii] = -area/M_PI*normal[ii]; - Ar1_2[ii] = Ar1[ii]/2.0; - for (jj=0; jj<3; jj++) { - nn[ii][jj] = normal[ii]*normal[jj]; - Ar2[ii][jj] = -area/M_PI*nn[ii][jj]; - Ar2_2[ii][jj] = Ar2[ii][jj]/2.0; - if (ii == jj) { - Ar3[ii][jj] = -area/M_PI*(2.0*nn[ii][jj] - 1.0); - } else { - Ar3[ii][jj] = -area/M_PI*(2.0*nn[ii][jj]); - } - Ar3_2[ii][jj] = Ar3[ii][jj]/2.0; - } - } -} - -// update Rho and S, generally for the case when they weren't originally input -void Facet::updateRhoS(double rhoSet, double sSet) -{ - rho = rhoSet; - - s = sSet; -} - -/*void Facet::getOpticalFourier(double* nnOut, double* Ar1Out, double* Ar2Out, double* Ar3Out, double* Ar1_2Out, double* Ar2_2Out, double* Ar3_2Out) -{ - nnOut = &nn[0][0]; - Ar1Out = &Ar1[0]; - Ar2Out = &Ar2[0][0]; - Ar3Out = &Ar3[0][0]; - Ar1_2Out = &Ar1_2[0]; - Ar2_2Out = &Ar2_2[0][0]; - Ar3_2Out = &Ar3_2[0][0]; -}*/ - -void Facet::getOpticalFourier(double nnOut[3][3], double Ar1Out[3], double Ar2Out[3][3], double Ar3Out[3][3], double Ar1_2Out[3], double Ar2_2Out[3][3], double Ar3_2Out[3][3]) -{ - int ii, jj; - - for (ii=0; ii<3; ii++) { - Ar1Out[ii] = Ar1[ii]; - Ar1_2Out[ii] = Ar1_2[ii]; - for (jj=0; jj<3; jj++) { - nnOut[ii][jj] = nn[ii][jj]; - Ar2Out[ii][jj] = Ar2[ii][jj]; - Ar3Out[ii][jj] = Ar3[ii][jj]; - Ar2_2Out[ii][jj] = Ar2_2[ii][jj]; - Ar3_2Out[ii][jj] = Ar3_2[ii][jj]; - } - } -} \ No newline at end of file diff --git a/lib/source/SRPModel.cpp b/lib/source/SRPModel.cpp deleted file mode 100644 index 1773c4d..0000000 --- a/lib/source/SRPModel.cpp +++ /dev/null @@ -1,868 +0,0 @@ -/** MIT License - -Copyright (c) 2014 Jay McMahon - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. -*/ - - -// -// SRPModel.cpp -// Spacecraft_SRP -// -// Created by Jay McMahon on 5/27/14. -// Copyright (c) 2014 Jay McMahon. All rights reserved. -// - -#include "SRPModel.h" -#include - -SRPModel::SRPModel(){ - return; -} - -SRPModel::SRPModel(double lambdaDel, double deltaDel, int MaxFourier, Body* bodyExist, int bounces, int numRefine){ - - // Assign body info - Spacecraft = *bodyExist; - - // Assign max order for Fourier coefficient - FourierOrder = MaxFourier; - - // Assign number of bounces - numBounces = bounces; - - // Assign how many refinement steps to take - numShadRefine = numRefine; - - // Set delta_s and lambda_s lists - int ii, jj; - int dsnum = 180/deltaDel + 1; - int lsnum = 360/lambdaDel + 1; - deltaSunList.resize(dsnum); - lambdaSunList.resize(lsnum); - riseLambda.resize(dsnum); - setLambda.resize(dsnum); - - for (ii=0; ii row(numf); - - for (ii=0; ii 0) { - Spacecraft.setViewRc(); - computeMulti(); - // Specular reflections - // if flag is set? or always - specularMulti(dsnum, lsnum); - } - - // Create coefficients array and fill - vector Frow(FourierOrder); - Coefficients.resize(dsnum); - shadowing.resize(lsnum); - shadowBoundsLow.resize(numf); - shadowBoundsHigh.resize(numf); - - for (ii=0; ii= 0){ - // Always illuminated - *rise = 2*M_PI; - *set = 0.0; - } - else{ - // Sometimes illuminated - nLong = atan2(normal[1],normal[0]); - arg = -(czo*czs)/(szo*szs); - dellon = acos(arg); - *rise = fmod((nLong + dellon), (2*M_PI)); - if (*rise < 0.0) *rise += 2*M_PI; - *set = fmod((nLong - dellon), (2*M_PI)); - if (*set < 0.0) *set += 2*M_PI; - } - return; -} - -void SRPModel::writeSRPCoeffs(int deltaWrite) -{ - int ii; - double* Ap; - double* Bp; - - // cout << deltaSunList[deltaWrite]*180.0/M_PI << " "; - Ap = Coefficients[deltaWrite][0].getA(); - // cout << *Ap << " " << *(Ap+1) << " " << *(Ap+2) << " "; - for (ii=1; ii* view; - vector adjRows; - vector adjCols; - vector < vector < double > > relVec; - vector < double > relVecR; - vector < vector < double > > relVecHat; - vector row(3); - double* rc1; - double* rc2; - double temp; - - // Get the count and required data for facets that see eachother - multiCount = 0; - for (ii=0; iisize(); jj++) { - if (view->at(jj) > ii) { - multiCount++; - } - } - } - adjRows.resize(multiCount); - adjCols.resize(multiCount); - relVec.resize(multiCount); - relVecR.resize(multiCount); - relVecHat.resize(multiCount); - int currentCount = 0; - for (ii=0; iisize(); jj++) { - if (view->at(jj) > ii) { - adjRows[currentCount] = ii; - adjCols[currentCount] = view->at(jj); - rc2 = Spacecraft.getFacetRc(view->at(jj)); - for (kk=0; kk<3; kk++) { - row[kk] = rc2[kk] - rc1[kk]; - } - relVec[currentCount] = row; - temp = sqrt(row[0]*row[0]+row[1]*row[1]+row[2]*row[2]); - relVecR[currentCount] = temp; - for (kk=0; kk<3; kk++) { - row[kk] /= temp; - } - relVecHat[currentCount] = row; - currentCount++; - } - } - } - - // Compute reflectivity matrix entries - vector reflMatRows(2*multiCount); - vector reflMatCols(2*multiCount); - vector reflMatVals(2*multiCount); - double theta_mid, deltaAng, pwrBounces; - for (ii=0; ii 4154){ - printf("i=%d\n",i); - printf("theta = %.5g\n",theta_mid); - printf("area col = %.5g\n",areas[adjCols[i]]); - printf("r = %.5g\n",relVecR[i]); - printf("rhat = %.5g %.5g %.5g\n",relVecHat[3*i],relVecHat[3*i+1],relVecHat[3*i+2]); - printf("norm col = %.5g %.5g %.5g\n",normals[3*adjCols[i]],normals[3*adjCols[i]+1],normals[3*adjCols[i]+2]); - printf("dAng = %.5g\n",deltaAng); - } - */ - - // Compute fraction from facet2 to facet1 - theta_mid = acos(-1.0 * dot(&relVecHat[ii][0], Spacecraft.getFacetNormal(adjCols[ii]))); - deltaAng = sqrt(Spacecraft.getFacetArea(adjRows[ii]) * dot(&relVecHat[ii][0], Spacecraft.getFacetNormal(adjRows[ii]))/(relVecR[ii]*relVecR[ii])); - reflMatRows[2*ii + 1] = adjCols[ii]; - reflMatCols[2*ii + 1] = adjRows[ii]; - reflMatVals[2*ii + 1] = (deltaAng/(2*M_PI))*sin(2.0*theta_mid)*sin(deltaAng); - - - // if only modeling 1 or 2 bounces, use the series power bounce to increase fidelity - if (numBounces < 3) { - // Compute power series for repeated reflections - pwrBounces = 1/(1-reflMatVals[2*ii]*reflMatVals[2*ii + 1]); - - // Put power series into reflectivity matrix - reflMatVals[2*ii] = pwrBounces*reflMatVals[2*ii]; - reflMatVals[2*ii + 1] = pwrBounces*reflMatVals[2*ii + 1]; - } - - } - - // printf("Refl[2 to 19] = %.5g\n",reflMatVals[0]); - // printf("Refl[19 to 2] = %.5g\n",reflMatVals[1]); - - // Compute all paths up to N bounces - // Currently, for the sake of getting it working, I'm only going to do 1 bounce - double tempvec[3], tempvec2[3]; - double Lambertian = 2.0/3.0; - double *tempn; - double *temprc; - - // initialize all entries with zeros - fMult.resize(numF); - tMult.resize(numF); - row[0] = 0.0; - row[1] = 0.0; - row[2] = 0.0; - for (ii=0; ii1) { - - // Make a copy of the single bounce fmult and tmult - vector < vector < double > > fMultLast(fMult); - vector < vector < double > > tMultLast(tMult); - - // Loop over number of bounces requested - int bounce = 2; - vector < vector < double > > fMultNew; - vector < vector < double > > tMultNew; - // initialize all entries with zeros - fMultNew.resize(numF); - tMultNew.resize(numF); - row[0] = 0.0; - row[1] = 0.0; - row[2] = 0.0; - for (ii=0; ii* view; - double uhat[3], vmag, pt[3]; - double hit_pnt[3], lightSource[3]; - double ftemp[3], ttemp[3], fi[3], mi[3]; - double uhat_new[3], scale_light; - double* rc1; - double *tempn; - int currF, currSpecCase; - double currDelS, currLamS; - Facet tempF, tempF2; - vector rowd; - vector < vector < double > > rowvd; - vector extraFacets; - - numF = Spacecraft.getNumFacets(); - - fSpec.resize(numF); - tSpec.resize(numF); - latSpec.resize(numF); - longSpec.resize(numF); - - for (ii=0; iiempty() && (tempF.getS()>0.0)) { - - // lat/long loop - for (mm=0; mm riseLambda[mm][ii]) && (currLamS < setLambda[mm][ii])) || ((riseLambda[mm][ii] > setLambda[mm][ii]) && ((currLamS > riseLambda[mm][ii]) || (currLamS < setLambda[mm][ii]))) ) { - continue; - } - lightSource[0] = cos(currDelS)*cos(currLamS); - lightSource[1] = cos(currDelS)*sin(currLamS); - lightSource[2] = sin(currDelS); - for (kk=0; kk<3; kk++) { - ftemp[kk] = 0.0; - ttemp[kk] = 0.0; - } - - bounce = 1; - scale_light = tempF.getRho() * tempF.getS(); - rc1 = Spacecraft.getFacetRc(currF); - while (bounce <= numBounces) { - // determine direction of ray to search - tempn = Spacecraft.getFacetNormal(currF); - vmag = dot(&lightSource[0], tempn); - // skip if dot product is too close to zero - neglible bouncing power and gives bad numerical results - if (vmag < 1e-15) { - test = -1; - } else { - for (kk=0; kk<3; kk++) { - uhat[kk] = -lightSource[kk] + 2*vmag*tempn[kk]; - } - vmag = sqrt(dot(&uhat[0], &uhat[0])); - for (kk=0; kk<3; kk++) { - uhat[kk] = uhat[kk]/vmag; - pt[kk] = rc1[kk] + 0.001*uhat[kk]; // move starting point 0.1% of the distance along the connecting line so the ray tracing check doesn't just return the starting facet - } - test = Spacecraft.ray_intersection(&pt[0], &uhat[0], &hit_pnt[0], &extraFacets); - } - // If there is a hit, and it hits from the front of the facet - if ((test != -1) && (dot(&uhat[0], Spacecraft.getFacetNormal(test))<-1e-15)) { - // Figure out the force/moment updates - // should just be a multiplier/factor I can plug into my coefficient calculations - // associate with each facet that receives a subsequent bounce - // computing f_i and m_i here; units of [area^2] - for (kk=0; kk<3; kk++) { - uhat_new[kk] = -uhat[kk]; - } - tempF = Spacecraft.getFacet(test); - compute_fi(&uhat_new[0], tempF, &fi[0]); - rc1 = Spacecraft.getFacetRc(test); - cross(rc1, &fi[0], &mi[0]); - for (kk=0; kk<3; kk++) { - ftemp[kk] += scale_light*fi[kk]; - ttemp[kk] += scale_light*mi[kk]; - } - - if (!extraFacets.empty()) { - // Go through any other facets hit at the same time - // at this point, we only compute the force from these extra facets, but bounce is only computed from the first hit facet - // I think what I want to do is an average of the slopes of the facets hit to compute the next bounce... also compute an average for updating lightSource? not sure I can do this with the dependence on currF though... - for (int qq=0; qq rowd; - int kk, lowhigh; - - // Compute shadowing information - for (int jj=0; jj shadowBoundsHigh[jj].size()) { - shadowBoundsHigh[jj].push_back(lambdaSunList[lsnum-1]); - } - } - } - - // Refine shadowing longitude boundaries if desired - if (numShadRefine > 0) { - int bound_count, shadFlag, mm, nn, testHit; - double newLow, newHigh, midp, maxDim; - double uhat[3], negu[3], pt[3]; - double* rc; - double hit_pnt[3]; - vector dummy; - maxDim = Spacecraft.getMaxDimVoxGrid(); - for (int jj=0; jj 0.0) { - shadFlag = 0; - newLow = shadowBoundsLow[jj][kk] - lambdaDel; - newHigh = shadowBoundsLow[jj][kk]; - for (mm=0; mmsize()>0) { - rc = Spacecraft.getFacetRc(jj); - for (nn=0; nn<3; nn++) { - pt[nn] = rc[nn] + 2.0*maxDim*uhat[nn]; - } - testHit = Spacecraft.ray_intersection(&pt[0], &negu[0], &hit_pnt[0], &dummy); - if (testHit != jj) { - // Facet was shadowed by another facet - if (shadFlag==0) { - newLow = midp; - } else if (shadFlag==1) { - newHigh = midp; - } - } else { - if (shadFlag==0) { - newHigh = midp; - } else if (shadFlag==1) { - newLow = midp; - } - } - } - } - if (shadFlag==0) { - shadowBoundsLow[jj][kk] = newHigh; - } else if (shadFlag==1) { - shadowBoundsHigh[jj][kk] = newLow; - } - } - if (shadowBoundsHigh[jj][kk] < 2.0*M_PI) { - shadFlag = 1; - newLow = shadowBoundsHigh[jj][kk]; - newHigh = shadowBoundsHigh[jj][kk] + lambdaDel; - for (mm=0; mmsize()>0) { - rc = Spacecraft.getFacetRc(jj); - for (nn=0; nn<3; nn++) { - pt[nn] = rc[nn] + 2.0*maxDim*uhat[nn]; - } - testHit = Spacecraft.ray_intersection(&pt[0], &negu[0], &hit_pnt[0], &dummy); - if (testHit != jj) { - // Facet was shadowed by another facet - if (shadFlag==0) { - newLow = midp; - } else if (shadFlag==1) { - newHigh = midp; - } - } else { - if (shadFlag==0) { - newHigh = midp; - } else if (shadFlag==1) { - newLow = midp; - } - } - } - } - if (shadFlag==0) { - shadowBoundsLow[jj][kk] = newHigh; - } else if (shadFlag==1) { - shadowBoundsHigh[jj][kk] = newLow; - } - } - } - } - } -} \ No newline at end of file diff --git a/lib/source/Shadow.cpp b/lib/source/Shadow.cpp deleted file mode 100644 index 15409fd..0000000 --- a/lib/source/Shadow.cpp +++ /dev/null @@ -1,168 +0,0 @@ -/** MIT License - -Copyright (c) 2014 Jay McMahon - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. -*/ - - -// -// Shadow.cpp -// Spacecraft_SRP -// -// Created by Jay McMahon on 5/27/14. -// Copyright (c) 2014 Jay McMahon. All rights reserved. -// - -#include "Shadow.h" -#include - -Shadow::Shadow() -{ - -} - -Shadow::Shadow(int numfacets, double delta, double lambda) -{ - deltaSun = delta; - lambdaSun = lambda; - - percentUnShadowed.resize(numfacets); -} - -void Shadow::ComputeShadowing(Body* Target, vector* riseLam, vector* setLam) -{ - // this is where we can use the inView information to simplify what facets to look through when testing to see if shadowed... will have to do comparison tests to see if it makes a significant difference (voxels might be fast enough) - - // to start, just ignore inView and test every case - int numf = Target->getNumFacets(); - int ii, jj; - int testHit; - vector dummy; - - // Define uhat - double uhat[3], negu[3], pt[3];//, rmag; - uhat[0] = cos(lambdaSun)*cos(deltaSun); - uhat[1] = sin(lambdaSun)*cos(deltaSun); - uhat[2] = sin(deltaSun); - for (jj=0; jj<3; jj++) { - negu[jj] = -1.0*uhat[jj]; - } - - double maxDim = Target->getMaxDimVoxGrid(); - - double* rc; -// double* rc2; - double hit_pnt[3]; -// double degBlock = M_PI/180.0; // would be better to input this on a per facet basis, giving a radius of a circle around rc that is within the circle, and then this can be computed to an angle based on facet-to-facet range. - - for (ii=0; iigetInView(ii)).size()>0) { - rc = Target->getFacetRc(ii); - // This logic starts "at the Sun" and looks toward the facet in quesiton - // Doesn't need to check for multiple facets as it should hit the center - for (jj=0; jj<3; jj++) { - pt[jj] = rc[jj] + 2.0*maxDim*uhat[jj]; - } - testHit = Target->ray_intersection(&pt[0], &negu[0], &hit_pnt[0], &dummy); - if (testHit != ii) { - // Facet was shadowed by another facet - percentUnShadowed[ii] = 0; - } - } - // Alternative is to use already computed inViewRc data; check each possible facet in view, and if the relative position is within X degrees of uhat, this facet is shadowed -// if ((*Target->getInViewRc(ii)).size()>0) { -// rc = Target->getFacetRc(ii); -// for (int kk=0; kk<(*Target->getInViewRc(ii)).size(); kk++ ){ -// rc2 = Target->getFacetRc((*Target->getInViewRc(ii))[kk]); -// for (jj=0; jj<3; jj++) { -// pt[jj] = rc2[jj] - rc[jj]; -// } -// rmag = sqrt(dot(pt,pt)); -// if (acos(dot(uhat, pt)/rmag) < degBlock) { -// // Facet was shadowed by another facet -// percentUnShadowed[ii] = 0; -// break; -// } -// } -// } - } else { - if (((*riseLam)[ii] < (*setLam)[ii]) && (lambdaSun > (*riseLam)[ii]) && (lambdaSun < (*setLam)[ii])) { - // Visible region wraps around 2PI, and Sun is outside of it so facet doesn't see Sun - percentUnShadowed[ii] = 0; - } else if (((*riseLam)[ii] > (*setLam)[ii]) && ((lambdaSun > (*riseLam)[ii]) || (lambdaSun < (*setLam)[ii]))) { - // Visible region is continuous, and Sun is outside of it so facet doesn't see Sun - percentUnShadowed[ii] = 0; - } else { - // Facet currently sees the sun - // Test to see if ray from Sun direction intersects this facet - // To add inView, do so here; start ray at the current facet and test only for intersection of the inView facets - if ((*Target->getInView(ii)).size()>0) { - rc = Target->getFacetRc(ii); - // This logic starts "at the Sun" and looks toward the facet in quesiton - // Doesn't need to check for multiple facets as it should hit the center - for (jj=0; jj<3; jj++) { - pt[jj] = rc[jj] + 2.0*maxDim*uhat[jj]; - } - testHit = Target->ray_intersection(&pt[0], &negu[0], &hit_pnt[0], &dummy); - if (testHit != ii) { - // Facet was shadowed by another facet - percentUnShadowed[ii] = 0; - } - } - // Alternative is to use already computed inViewRc data; check each possible facet in view, and if the relative position is within X degrees of uhat, this facet is shadowed -// if ((*Target->getInViewRc(ii)).size()>0) { -// rc = Target->getFacetRc(ii); -// for (int kk=0; kk<(*Target->getInViewRc(ii)).size(); kk++ ){ -// rc2 = Target->getFacetRc((*Target->getInViewRc(ii))[kk]); -// for (jj=0; jj<3; jj++) { -// pt[jj] = rc2[jj] - rc[jj]; -// } -// rmag = sqrt(dot(pt,pt)); -// if (acos(dot(uhat, pt)/rmag) < degBlock) { -// // Facet was shadowed by another facet -// percentUnShadowed[ii] = 0; -// break; -// } -// } -// } - } - } - } -} - -int Shadow::getPcntLit(int fnum) -{ - return percentUnShadowed[fnum]; -} - -double Shadow::dot ( double* a, double* b) -{ - return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; - -} \ No newline at end of file diff --git a/lib/source/VoxelGrid.cpp b/lib/source/VoxelGrid.cpp deleted file mode 100644 index 90fb5a6..0000000 --- a/lib/source/VoxelGrid.cpp +++ /dev/null @@ -1,783 +0,0 @@ -/** MIT License - -Copyright (c) 2014 Jay McMahon - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. -*/ - - -// -// VoxelGrid.cpp -// Spacecraft_SRP -// -// Created by Jay McMahon on 5/19/14. -// Copyright (c) 2014 Jay McMahon. All rights reserved. -// - -#include "VoxelGrid.h" -#include -#include -#include - -using namespace std; - -VoxelGrid::VoxelGrid() -{ - VoxelGrid(1.0,1.0,1.0,1); -} - -VoxelGrid::VoxelGrid(double xmax_in, double ymax_in, double zmax_in, int N_in){ - N = N_in; - dx = 2.0*xmax_in/(double)N; - dy = 2.0*ymax_in/(double)N; - dz = 2.0*zmax_in/(double)N; - - // Resize vectors in class - int N3 = N*N*N; - int N2 = N*N; - xmax.resize(N3); - xmin.resize(N3); - ymax.resize(N3); - ymin.resize(N3); - zmax.resize(N3); - zmin.resize(N3); - posx.resize(N3); - negx.resize(N3); - posy.resize(N3); - negy.resize(N3); - posz.resize(N3); - negz.resize(N3); - - flist.resize(N3); - vlist.resize(N3); - fnum.resize(N3); - vnum.resize(N3); - - vector row; - - - int vox_ind = 0; - - int ii, jj, kk; - - for (ii = 0; ii* facets, vector* vertices, int numVerts, vector< vector >* f2vlist, vector* fdata) -{ - int ii; - vector voxel_list; - vector::iterator iter_vl; - vector< vector > v2voxlist; // local list of voxels each vertex is in - v2voxlist.resize(numVerts); - - //num_vertices = size(vertices,1); - - //f_2_v_list = facets_per_vertex(facets,num_vertices); - ii = 0; - for (vector::iterator iter = vertices->begin(); iter != vertices->end(); iter+=3) { - - double v1[3]; - v1[0] = *iter; - v1[1] = *(iter+1); - v1[2] = *(iter+2); - - find_pt_in_voxel(&v1[0], &voxel_list); - v2voxlist[ii] = voxel_list; - - for (iter_vl = voxel_list.begin(); iter_vl != voxel_list.end(); ++iter_vl) { - - // add vertex and associated facets to this voxel's list - vlist[*iter_vl].push_back(ii); - vnum[*iter_vl]++; -// flist[*iter_vl].reserve(flist[*iter_vl].size() + f2vlist[ii].size()); - copy((*f2vlist)[ii].begin(), (*f2vlist)[ii].end(), back_inserter(flist[*iter_vl])); - sort(flist[*iter_vl].begin(),flist[*iter_vl].end()); - flist[*iter_vl].erase(unique(flist[*iter_vl].begin(),flist[*iter_vl].end()),flist[*iter_vl].end()); - fnum[*iter_vl] = flist[*iter_vl].size(); - } - voxel_list.erase(voxel_list.begin(),voxel_list.end()); - ii++; - } - - // Check each facet's vertices to see if the facet spans more than one voxel - vector facet_voxs, checkList; - int temp = 0, temp2 = 0, temp3 = 0, jj, kk, mm; - int N2 = N*N; - int N3 = N2*N; - int fxmin,fxmax,fymin,fymax,fzmin,fzmax; - ii=0; // index of facets (for flist/fnum) - for (vector::iterator iterf = facets->begin(); iterf != facets->end(); iterf+=3) { - // Recall that facetList contains the numbers from the input file, so the value referenced by iterf needs to have 1 subtracted to start at index 0 - facet_voxs = v2voxlist[*(iterf)-1]; - copy(v2voxlist[*(iterf+1)-1].begin(), v2voxlist[*(iterf+1)-1].end(), back_inserter(facet_voxs)); - copy(v2voxlist[*(iterf+2)-1].begin(), v2voxlist[*(iterf+2)-1].end(), back_inserter(facet_voxs)); - sort(facet_voxs.begin(), facet_voxs.end()); - facet_voxs.erase(unique(facet_voxs.begin(), facet_voxs.end()), facet_voxs.end()); - if (facet_voxs.size()>1) { - if (facet_voxs.size()==2) { - // if there are only two voxels and they share a face that is all possible - temp = abs(facet_voxs[1] - facet_voxs[0]); - if ((temp==1) || (temp==N) || (temp == N2)) { - ii++; - continue; - } - } - // otherwise, check all voxels in between the max/min on each axis to see if the triangle intersects - fxmin = N3; - fxmax = -1; - fymin = N3; - fymax = -1; - fzmin = N3; - fzmax = -1; - for (jj=0; jj::iterator iterc = checkList.begin(); iterc != checkList.end(); ++iterc) { - // function returns true if an intersection - if (intersectTriVox(*iterc, &(*vertices)[(*fdata)[ii].getVert(0)], &(*vertices)[(*fdata)[ii].getVert(1)], &(*vertices)[(*fdata)[ii].getVert(2)], (*fdata)[ii].getNormal())) { - // there is an intersection; add facet to appropriate voxel - flist[*iterc].push_back(ii); - fnum[*iterc]++; - } - } - } - //facet_voxs.erase(facet_voxs.begin(),facet_voxs.end()); - ii++; - } - - // Remove duplicate facet entries - int vox_ind = 0; - for (ii = 0; ii* voxel_list) -{ - // Test to ensure the pt is inside the total voxel_grid box - if ((pt[0] < xmin[0]) || (pt[0] > abs(xmin[0])) || (pt[1] < ymin[0]) || (pt[1] > abs(ymin[0])) || (pt[2] < zmin[0]) || (pt[2] > abs(zmin[0]))) { - return; - } - - // Find voxels this pt is inside of (or on the border) - // I know that the first voxel is at the negative end of all three axes - // Account for numerical rounding by comparing to floor - double num_accuracy = 1e-14; - double xgrid, ygrid, zgrid; - int xind, yind, zind; - - xgrid = (pt[0] - xmin[0])/dx; - if (abs(xgrid - round(xgrid)) < num_accuracy) { - xgrid = round(xgrid); - } - xind = (int)floor(xgrid); - if (xind == N) { - xind--; - } - ygrid = (pt[1] - ymin[0])/dy; - if (abs(ygrid - round(ygrid)) < num_accuracy) { - ygrid = round(ygrid); - } - yind = (int)floor(ygrid); - if (yind == N) { - yind--; - } - zgrid = (pt[2] - zmin[0])/dz; - if (abs(zgrid - round(zgrid)) < num_accuracy) { - zgrid = round(zgrid); - } - zind = (int)floor(zgrid); - if (zind == N) { - zind--; - } - - voxel_list->push_back(xind + yind*N + zind*N*N); - - // check neighboring voxels - if ((xgrid == xind) && (xind != 0)) { - voxel_list->push_back((*voxel_list)[0] - 1); - if ((ygrid == yind) && (yind != 0)) { - voxel_list->push_back((*voxel_list)[0] - N - 1); - if ((zgrid == zind) && (zind != 0)) { - voxel_list->push_back((*voxel_list)[0] - N*N - N - 1); - } - } - if ((zgrid == zind) && (zind != 0)) { - voxel_list->push_back((*voxel_list)[0] - N*N - 1); - } - } - if ((ygrid == yind) && (yind != 0)) { - voxel_list->push_back((*voxel_list)[0] - N); - if ((zgrid == zind) && (zind != 0)) { - voxel_list->push_back((*voxel_list)[0] - N*N - N); - } - } - if ((zgrid == zind) && (zind != 0)) { - voxel_list->push_back((*voxel_list)[0] - N*N); - } - - return; -} - -int VoxelGrid::ray_intersect_voxel(double* pt, double* uhat, vector* facets, vector* vertices, vector * facetData, double* hit_pnt, vector* extraFacets) -{ - int hit_facet = -1; - vector voxel_list; - double tX, tY, tZ, t; - int current_voxel; - - // Find which voxel ray starts in - find_pt_in_voxel(pt, &voxel_list); - - // If starts outside voxel_grid, find where the ray intersects - while (voxel_list.size() == 0) { - // Figure out where first intersection occurs - t = 0.0; - if (abs(uhat[0]) < 1e-16) { - tX = -1.0; - } else { - if (pt[0] > 0) { - tX = (abs(xmin[0]) - pt[0])/uhat[0]; - } else { - tX = (xmin[0] - pt[0])/uhat[0]; - } - } - if (tX > 0.0) { - t = tX; - } - if (abs(uhat[1]) < 1e-16) { - tY = -1.0; - } else { - if (pt[1] > 0.0) { - tY = (abs(ymin[0]) - pt[1])/uhat[1]; - } else { - tY = (ymin[0] - pt[1])/uhat[1]; - } - } - if ((tY > 0.0) && ((t==0.0) || (tY < t))) { - t = tY; - } - if (abs(uhat[2]) < 1e-16) { - tZ = -1.0; - } else { - if (pt[2] > 0.0) { - tZ = (abs(zmin[0]) - pt[2])/uhat[2]; - } else { - tZ = (zmin[0] - pt[2])/uhat[2]; - } - } - if ((tZ > 0.0) && ((t==0.0) || (tZ < t))) { - t = tZ; - } - - // If ray is pointed away from the voxel_grid, return no intersection - if (t == 0.0) { - return hit_facet; // returns -1 if no intersection - } - - // Update pt - pt[0] = pt[0] + t*uhat[0]; - pt[1] = pt[1] + t*uhat[1]; - pt[2] = pt[2] + t*uhat[2]; - - // Check to see if this intersect a voxel now - find_pt_in_voxel(pt, &voxel_list); - - } - - if (voxel_list.size() > 1) { - // the ray starts on a border - figure out which of the list is the - // voxel the ray is going to travel through - // It is possible that the ray will stay on the border of two voxels - - // in that case, just choose one as the equalities for testing facets - // will work it out appropriately. - double stepsize = min(min(dx, dy), dz)/10.0; - double pt2[3]; - vector voxel_list2; - pt2[0] = pt[0] + stepsize*uhat[0]; - pt2[1] = pt[1] + stepsize*uhat[1]; - pt2[2] = pt[2] + stepsize*uhat[2]; - find_pt_in_voxel(pt2, &voxel_list2); - - if ((abs(uhat[0])< 1e-13) || (abs(uhat[1])< 1e-13) || (abs(uhat[2])< 1e-13)) { - // ray will travel on some border - decide appropriate facet to test - while (find(voxel_list.begin(),voxel_list.end(),voxel_list2[0])==voxel_list.end()) { - stepsize = stepsize/10.0; - pt2[0] = pt[0] + stepsize*uhat[0]; - pt2[1] = pt[1] + stepsize*uhat[1]; - pt2[2] = pt[2] + stepsize*uhat[2]; - voxel_list2.erase(voxel_list2.begin(),voxel_list2.end()); - find_pt_in_voxel(pt2, &voxel_list2); - } - current_voxel = voxel_list2[0]; - - } else { - // ray will enter one of the current voxels uniquely - find it - // ensure that we find one unique voxel, and that it is one of the originally bordered voxels - while ((voxel_list2.size() > 1) || (find(voxel_list.begin(),voxel_list.end(),voxel_list2[0])==voxel_list.end())) { - stepsize = stepsize/10.0; - pt2[0] = pt[0] + stepsize*uhat[0]; - pt2[1] = pt[1] + stepsize*uhat[1]; - pt2[2] = pt[2] + stepsize*uhat[2]; - voxel_list2.erase(voxel_list2.begin(),voxel_list2.end()); - find_pt_in_voxel(pt2, &voxel_list2); - } - current_voxel = voxel_list2[0]; - } - - } else { - current_voxel = voxel_list[0]; - } - - // Get step sizes to cross current voxel and full voxels in each direction - double tDeltaX, tDeltaY, tDeltaZ; - double tMaxX, tMaxY, tMaxZ; - - if (abs(uhat[0]) < 1e-16) { - tDeltaX = 2.0*dx*N; - tMaxX = 2.0*dx*N; - } else { - tDeltaX = dx/abs(uhat[0]); - if (uhat[0] > 0.0) { - tMaxX = (xmax[current_voxel] - pt[0])/uhat[0]; - } else { - tMaxX = (xmin[current_voxel] - pt[0])/uhat[0]; - } - } - if (abs(uhat[1]) < 1e-16) { - tDeltaY = 2.0*dy*N; - tMaxY = 2.0*dy*N; - } else { - tDeltaY = dy/abs(uhat[1]); - if (uhat[1] > 0.0) { - tMaxY = (ymax[current_voxel] - pt[1])/uhat[1]; - } else { - tMaxY = (ymin[current_voxel] - pt[1])/uhat[1]; - } - } - if (abs(uhat[2]) < 1e-16) { - tDeltaZ = 2.0*dz*N; - tMaxZ = 2.0*dz*N; - } else { - tDeltaZ = dz/abs(uhat[2]); - if (uhat[2] > 0.0) { - tMaxZ = (zmax[current_voxel] - pt[2])/uhat[2]; - } else { - tMaxZ = (zmin[current_voxel] - pt[2])/uhat[2]; - } - } - - // Traverse through the voxel grid - int N3 = N*N*N; - while (current_voxel != N3) { - - // Conduct test - if the test is satisfied then return - // Note that we have to make sure the test function returns an empty - // hit_facet if there is no hit - // If current voxel doesn't contain any facets, just traverse to next - // voxel. No need to take the time to call the function - if (fnum[current_voxel] != 0) { - hit_facet = test_intersection(pt, uhat, facets, vertices, facetData, hit_pnt, current_voxel, extraFacets); - - if (hit_facet != -1) { - return hit_facet; - } - } - - // Traverse to next voxel - if (tMaxX < tMaxY) { - if (tMaxX < tMaxZ) { - if (uhat[0] > 0.0) { - current_voxel = posx[current_voxel]; - } else { - current_voxel = negx[current_voxel]; - } - tMaxX= tMaxX + tDeltaX; - } else { - if (uhat[2] > 0.0) { - current_voxel = posz[current_voxel]; - } else { - current_voxel = negz[current_voxel]; - } - tMaxZ= tMaxZ + tDeltaZ; - } - } else { - if (tMaxY < tMaxZ) { - if (uhat[1] > 0.0) { - current_voxel = posy[current_voxel]; - } else { - current_voxel = negy[current_voxel]; - } - tMaxY= tMaxY + tDeltaY; - } else { - if (uhat[2] > 0.0) { - current_voxel = posz[current_voxel]; - } else { - current_voxel = negz[current_voxel]; - } - tMaxZ= tMaxZ + tDeltaZ; - } - } - -// if (hit_facet != -1) { -// current_voxel = N3; -// } - } - - // If we made it here, it is because the ray left the grid without intersecting any facets - // hit_facet will still be = -1, and hit_pnt will be empty - return hit_facet; -} - -int VoxelGrid::test_intersection(double* pt, double* uhat, vector* facets, vector* vertices, vector * facetData, double* hit_pnt, int current_voxel, vector* extraFacets) -{ - int hit_facet = -1; - - // Just need to find the first facet hit, so have to look at all of the - // facets in the group, and select the shortest intersection if there are - // more than one. - - // Pay no attention to which "side" of the triangle is intersected; just - // pick the first one along the path - - // Note that I could change this to output t, which is the range to the - // intersection along the ray - - int num_facets = fnum[current_voxel]; - int ii, current_facet; - - double tmin = 10.0*N*max(max(dx,dy),dz); // larger than the total voxel grid to start - double pt_int[3], t, temp1, temp2; - double* nhat; - double* v1; - double* v2; - double* v3; - bool hit; - - for (ii = 0; ii < num_facets; ii++) { - - // Solve for t in ray equation that gives a point in the triangle plane - // dot((pt + t.*uhat) - v, normal) == 0 where v is one of the vertices in the triangle - // dot(pt - v,normal) + t*dot(uhat,normal) == 0 - // t = dot(v - pt,normal)/dot(uhat,normal) - current_facet = flist[current_voxel][ii]; - nhat = (*facetData)[current_facet].getNormal(); - v1 = &(*vertices)[(*facetData)[current_facet].getVert(0)]; - v2 = &(*vertices)[(*facetData)[current_facet].getVert(1)]; - v3 = &(*vertices)[(*facetData)[current_facet].getVert(2)]; - pt_int[0] = v1[0] - pt[0]; - pt_int[1] = v1[1] - pt[1]; - pt_int[2] = v1[2] - pt[2]; - temp1 = dot(&pt_int[0],nhat); - temp2 = dot(uhat,nhat); - t = temp1/temp2; - - // Ray is directed; don't travel backwards to an intersection - if (t >= 0.0) { - pt_int[0] = pt[0] + t*uhat[0]; - pt_int[1] = pt[1] + t*uhat[1]; - pt_int[2] = pt[2] + t*uhat[2]; - - // test this point to see if its inside the current facet - hit = in_triangle_3D(v1, v2, v3, &pt_int[0], nhat); - - if (hit) { - if (t < tmin) { - hit_facet = current_facet; - hit_pnt[0] = pt_int[0]; - hit_pnt[1] = pt_int[1]; - hit_pnt[2] = pt_int[2]; - tmin = t; - } - else if (abs(t - tmin) < 1e-14) { - extraFacets->push_back(current_facet); - } - } - } - } - - return hit_facet; -} - -void VoxelGrid::inverse ( double* answer, double* x, double* y, double* n) -{ - double denom; - - denom = x[0]*(y[1]*n[2]-n[1]*y[2]) - y[0]*(x[1]*n[2]-n[1]*x[2]) + n[0]*(x[1]*y[2]-y[1]*x[2]); - - answer[0] = (y[1]*n[2]-n[1]*y[2])/denom; - answer[1] = (y[2]*n[0]-n[2]*y[0])/denom; - answer[2] = (y[0]*n[1]-n[0]*y[1])/denom; - answer[3] = (x[2]*n[1]-n[2]*x[1])/denom; - answer[4] = (x[0]*n[2]-n[0]*x[2])/denom; - answer[5] = (x[1]*n[0]-n[1]*x[0])/denom; - // answer[6] = (y[2]*x[1]-x[2]*y[1])/denom; - // answer[7] = (y[0]*x[2]-x[0]*y[2])/denom; - // answer[8] = (y[1]*x[0]-x[1]*y[0])/denom; - - return; -} - -bool VoxelGrid::in_triangle_3D( double *A, double *B, double *C, double *P, double *nhat) -{ - // In triangle stuff - double v0[3], v1[3], v2[3], answer[6], uv[2]; - int i; - - for (i=0; i<3; i++){ - v0[i] = C[i] - A[i]; - v1[i] = B[i] - A[i]; - v2[i] = P[i] - A[i]; - } - - inverse( &answer[0], &v1[0], &v0[0], nhat ); - - uv[0] = answer[0]*v2[0] + answer[1]*v2[1] + answer[2]*v2[2]; - uv[1] = answer[3]*v2[0] + answer[4]*v2[1] + answer[5]*v2[2]; - - return ((uv[0] >= -1e-14) && (uv[1] >= -1e-14) && (uv[0] + uv[1] <= (1.0+1e-14))); -} - -double VoxelGrid::dot ( double* a, double* b) -{ - return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; - -} - -int VoxelGrid::getNumFinVox(int VoxNum) -{ - return fnum[VoxNum]; -} - -int VoxelGrid::getNumVinVox(int VoxNum) -{ - return vnum[VoxNum]; -} - -double VoxelGrid::getMaxDim() -{ - double x = 0.5*dx*N; - double y = 0.5*dy*N; - double z = 0.5*dz*N; - - return sqrt(x*x + y*y + z*z); -} - -bool VoxelGrid::intersectTriVox(int voxNum, double* v1, double* v2, double* v3, double* nhat) -{ - // test for an intersection of the current voxel and facet using the separating axis theorem - // returns true if there IS an intersection, false if no intersection - double axis[3], rtemp; - double hx = dx/2.0; - double hy = dy/2.0; - double hz = dz/2.0; - double center[3], vc1[3], vc2[3], vc3[3]; - center[0] = xmin[voxNum] + hx; - center[1] = ymin[voxNum] + hy; - center[2] = zmin[voxNum] + hz; - int ii; - // translate facet to voxel centered frame - for (ii = 0; ii<3; ii++) { - vc1[ii] = v1[ii] - center[ii]; - vc2[ii] = v2[ii] - center[ii]; - vc3[ii] = v3[ii] - center[ii]; - } - - // First test the three voxel face normal directions, which are assumed to be axis aligned - axis[0] = 1.0; axis[1] = 0.0; axis[2] = 0.0; - if (satAABBTri(hx, &axis[0], &vc1[0], &vc2[0], &vc3[0])) { - return false; - } - axis[0] = 0.0; axis[1] = 1.0; axis[2] = 0.0; - if (satAABBTri(hy, &axis[0], &vc1[0], &vc2[0], &vc3[0])) { - return false; - } - axis[0] = 0.0; axis[1] = 0.0; axis[2] = 1.0; - if (satAABBTri(hz, &axis[0], &vc1[0], &vc2[0], &vc3[0])) { - return false; - } - - // Test facet normal - rtemp = hx*abs(nhat[0]) + hy*abs(nhat[1]) + hz*abs(nhat[2]); - if (satAABBTri(rtemp, nhat, &vc1[0], &vc2[0], &vc3[0])) { - return false; - } - - // Test facet edges crossed with voxel edges - double f0[3], f1[3], f2[3], e[3]; - for (ii = 0; ii<3; ii++) { - f0[ii] = v2[ii]-v1[ii]; - f1[ii] = v3[ii]-v2[ii]; - f2[ii] = v1[ii]-v3[ii]; - } - - e[0]= 0.0; e[1] = 0.0; e[2] = 0.0; - for (int jj=0; jj<3; jj++) { - e[jj] = 1.0; - cross(&e[0], &f0[0], &axis[0]); - rtemp = hx*abs(axis[0]) + hy*abs(axis[1]) + hz*abs(axis[2]); - if (satAABBTri(rtemp, &axis[0], &vc1[0], &vc2[0], &vc3[0])) { - return false; - } - e[jj] = 0.0; - } - for (int jj=0; jj<3; jj++) { - e[jj] = 1.0; - cross(&e[0], &f1[0], &axis[0]); - rtemp = hx*abs(axis[0]) + hy*abs(axis[1]) + hz*abs(axis[2]); - if (satAABBTri(rtemp, &axis[0], &vc1[0], &vc2[0], &vc3[0])) { - return false; - } - e[jj] = 0.0; - } - for (int jj=0; jj<3; jj++) { - e[jj] = 1.0; - cross(&e[0], &f2[0], &axis[0]); - rtemp = hx*abs(axis[0]) + hy*abs(axis[1]) + hz*abs(axis[2]); - if (satAABBTri(rtemp, &axis[0], &vc1[0], &vc2[0], &vc3[0])) { - return false; - } - e[jj] = 0.0; - } - - // No separating axis found, so they intersect - return true; - -} - -bool VoxelGrid::satAABBTri(double rBox, double* axis, double* v1, double* v2, double* v3) -{ - // returns TRUE if there IS a separating axis, indicating NO INTERSECTION - - double p1, p2, p3; - p1 = dot(axis, v1); - p2 = dot(axis, v2); - p3 = dot(axis, v3); - - if ((min(p1,min(p2,p3)) > rBox) || (max(p1,max(p2,p3)) < -rBox)) { - return true; - } else { - return false; - } -} - -void VoxelGrid::cross(double* a, double* b, double* c) -{ - c[0] = a[1]*b[2] - a[2]*b[1]; - c[1] = a[2]*b[0] - a[0]*b[2]; - c[2] = a[0]*b[1] - a[1]*b[0]; - return; -} \ No newline at end of file