Моделирование воздействия теплового излучения на элементы космического аппарата

vertex3(), v3d_dist() {

p3d_centr.x = 0;

p3d_centr.y = 0;

p3d_centr.z = 0;

}

//coords of objects

Vector3d v3d_normal;

Point3d vertex1;

Point3d vertex2;

Point3d vertex3;

Point3d p3d_centr;

Vector3d v3d_dist;

// area of a triangle

//float area;

float Area();

//triangles must know their radiosity

float sun_direct_B = 0;

float sun_reflected_B = 0;

float moon_temperature_B = 0;

float selfreflection_B = 0;

//albedo

float p = 0;

//dot product of normals of 2 triangle objects to define their visability

float operator * (const Model_Triangle &v) const {

return (v3d_normal.x * v.v3d_normal.x +

v3d_normal.y * v.v3d_normal.y +

v3d_normal.z * v.v3d_normal.z);

}

Model_Triangle& operator=(Model_Triangle& o){ //overload '=' to set

Triangles

if (this != &o){

v3d_normal = o.v3d_normal;

vertex1 = o.vertex1;

vertex2 = o.vertex2;

vertex3 = o.vertex3;

}

return *this;

}

~Model_Triangle();

};

Приложение 4

make_scene.h

#pragma once

#include "ModelVertext.h"

#include <deque>

#include <iostream>

#include <fstream>

#include <iomanip>

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <vector>

#include <deque>

#include <string>

#include <cstring>

using namespace std;

//make SCENE

void make_scene(deque<Model_Triangle>& SCENE, fstream& fin, float p);

float char_to_float(string& char_buffer,

int& n);

//compute number triangles

int calc_N(string *ref);

func.h

#pragma once

#include "assert.H"

#include <deque>

#include <iostream>

#include <fstream>

#include <iomanip>

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <vector>

#include <deque>

#include <string>

#include <cstring>

#include <algorithm>

#include "TRIANGLES.H"

#include "POINT.H"

#include "TRANSFORMATION.h"

//cspice

#include "C:\naif\cspice32\include\SpiceUsr.h"

using namespace std;

//get a vector of distance between centrs of 2 triangles

Vector3d Distance(Model_Triangle& _i, Model_Triangle& _j);

//Norm a vector

float Norm_v3d(Vector3d& _v3d);

//return dot product of 2 vectors

float dot_product(Vector3d& a, Vector3d& b);

float dot_product(Vector3d& a, Point3d& b);

//return alpha vetween two triangles

float compute_SURF_alpha(Vector3d&,

Vector3d&,

Vector3d _ij_distance,

float(*dot_product)(Vector3d&, Vector3d&),

float(*Norm_v3d)(Vector3d&),

float _j_Area);

//read matrix from file and solve with G-Z method

void zeidel(deque<Model_Triangle>& SCENE, vector< vector < float > >&

reflection_matrix,

vector<float>& x,

int columns,

float(*char_to_float1)(string&),

bool(*converge)(vector<float>& xk, vector<float>& xkp, int n));

//define convergens of G-Z

bool converge(vector<float>& xk, vector<float>& xkp, int n);

//when read a martix realease a char buffer and convert to float

float char_to_float1(string& char_buffer);

void make_matrix(deque<Model_Triangle>& SCENE, vector< vector <

float > > &reflection_matrix);

int find_min(vector<float> t);

void cmp_fluxes_on_surf(deque<Model_Triangle>& SURF, SpiceDouble

et);

void cmp_model_flux_from_surf(deque<Model_Triangle>& SCENE,

deque<Model_Triangle>& SURF);

void cmp_direct_flux_on_box(deque<Model_Triangle>& BOX,

SpiceDouble et);

void cmp_direct_flux_on_model(deque<Model_Triangle>& SCENE,

SpiceDouble et);

bool is_oculted_moon(SpiceDouble et);

void find_sun();

Приложение 5

transformation.h

#pragma once

#include "TRIANGLES.h"

#include "FUNC.h"

#include <vector>

#include <iostream>

#include <cmath>

void ortho_project(Model_Triangle& triangle);

void translate(const Point3d& camera_pos, Model_Triangle& triangle);

Vector3d normalize(Vector3d vec);

void rotate(const Point3d& _pos, Model_Triangle& triangle, float angle);

bool isintersect(Model_Triangle& target, Vector3d& D, Vector3d& dist, int

n);

bool isintersectsun(Model_Triangle& target, Vector3d& D,

std::vector<float>& dist, int n);

void cam_view(const Point3d& _camera_pos,

const Vector3d _camera_normal,

Model_Triangle& triangle);

Приложение 6

triangles.cpp

#include "assert.H"

#include <iostream>

#include <vector>

#include "TRIANGLES.H"

using namespace std;

//initilize a triangle object

Model_Triangle::Model_Triangle(std::vector< std::vector<float> >& v_)

: v3d_normal(v_[0]), vertex1(v_[1]), vertex2(v_[2]), vertex3(v_[3]),

v3d_dist() {

//itillize barycentr

p3d_centr.x = (vertex1.x + vertex2.x + vertex3.x) / 3;

p3d_centr.y = (vertex1.y + vertex2.y + vertex3.y) / 3;

p3d_centr.z = (vertex1.z + vertex2.z + vertex3.z) / 3;

}

//compute area method

float Model_Triangle::Area() {

float Sx = ((vertex2.y - vertex1.y)*

(vertex3.z - vertex1.z) -

(vertex3.y - vertex1.y)*

(vertex2.z - vertex1.z)) / 2;

float Sy = ((vertex1.x * vertex3.z +

vertex2.z * vertex3.x +

vertex2.x * vertex1.z) -

(vertex1.x * vertex2.z +

vertex2.x * vertex3.z +

vertex3.x * vertex1.z)) / 2;

float Sz = ((vertex1.x * vertex2.y +

vertex1.y * vertex3.x +

vertex2.x * vertex3.y) -

(vertex3.x * vertex2.y +

vertex3.y * vertex1.x +

vertex2.x * vertex1.y)) / 2;

float area = sqrt(pow(Sx, 2.0) + pow(Sy, 2.0) + pow(Sz, 2.0));

return area;

}

Model_Triangle::~Model_Triangle() {

}

Приложение 7

make_scene.cpp

#define _CRT_SECURE_NO_WARNINGS

#include <iostream>

#include <fstream>

#include <iomanip>

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <vector>

#include <deque>

#include <string>

#include <cstring>

#include "TRIANGLES.H"

#include "POINT.H"

using namespace std;

int global_counter;

//this function releases the char buffer and make it content an element of an

array of floats

float char_to_float(string& char_buffer,

int& n);

//this function reads symbols from .STL, make an object of class triangle

//and then put an object into the que

void make_scene(deque<Model_Triangle>& SCENE, fstream& fin, float p)

{

global_counter = 0;

int tr_it0 = 0;

int tr_it1 = 0;

int tr_it2 = 0;

string::iterator i;

float triangle[4][3];

string char_buffer;

int n = 0;

char c;

//main loop which reads file symbol by symbol until eof()

while (!fin.eof()) {

c = fin.get();

if (c == char(32)) { //check end of line and push array if true

triangle[tr_it1][tr_it2] = char_to_float(char_buffer, n);

if (fabs(triangle[tr_it1][tr_it2]) < 0.00001) {

triangle[tr_it1][tr_it2] = 0;

}

char_buffer.erase(char_buffer.begin(), char_buffer.end());

tr_it2++;

n++;

continue;

}

if (c == char(95)){ //check end of a word (coord) and if true push array

triangle[tr_it1][tr_it2] = char_to_float(char_buffer, n);

if (fabs(triangle[tr_it1][tr_it2]) < 0.00001) {

triangle[tr_it1][tr_it2] = 0;

}

char_buffer.erase(char_buffer.begin(), char_buffer.end());

fin.ignore(16, char(10));

n++;

tr_it2 = 0;

tr_it1++;

if (tr_it1 == 4) {

tr_it1 = 0;

//make a [4][3] stl vector of verteces coordinates and normal

vector< vector < float > > tr_vec;

for (int i = 0; i < 4; i++){

vector<float> tr_vec_pt;

tr_vec.push_back(tr_vec_pt);

for (int j = 0; j < 3; j++){

tr_vec[i].push_back(triangle[i][j]);

}

}

//place a triangle object in que

SCENE.emplace_back(tr_vec);

SCENE[global_counter].p = p; //set albedo

global_counter++;

}

continue;

}

char_buffer.push_back(c);

}

}

float char_to_float(string& char_buffer,

int& n) {

float fl = 0;

char* temp = new char[char_buffer.length() + 1];

std::strcpy(temp, char_buffer.c_str());

fl = atof(temp);

return fl;

delete[] temp;

}

Приложение 8

transformation.cpp

#include "TRANSFORMATION.h"

#include <vector>

#include <iostream>

using namespace std;

Vector3d normalize(Vector3d vec) {

float length = sqrt(pow(vec.x, 2.0) + pow(vec.y, 2.0) + pow(vec.z, 2.0));

vec.x = vec.x / length;

vec.y = vec.y / length;

vec.z = vec.z / length;

return vec;

}

void ortho_project(Model_Triangle& triangle) {

float width = 524;

float height = 306;

float Zfar = 300;

float Znear = 0;

triangle.vertex1.x *= 2 / width;

triangle.vertex1.y *= 2 / height;

triangle.vertex1.z *= (-(2 / (Zfar - Znear)) - ((Zfar + Znear) / (Zfar -

Znear)));

triangle.vertex2.x *= 2 / width;

triangle.vertex2.y *= 2 / height;

triangle.vertex2.z *= (-(2 / (Zfar - Znear)) - ((Zfar + Znear) / (Zfar -

Znear)));

triangle.vertex3.x *= 2 / width;

triangle.vertex3.y *= 2 / height;

triangle.vertex3.z *= (-(2 / (Zfar - Znear)) - ((Zfar + Znear) / (Zfar -

Znear)));

triangle.v3d_normal.x *= 2 / width;

triangle.v3d_normal.y *= 2 / height;

triangle.v3d_normal.z *= (-(2 / (Zfar - Znear)) - ((Zfar + Znear) / (Zfar -

Znear)));

}

//translation

void translate(const Point3d& new_pos, Model_Triangle& triangle) {

triangle.vertex1.x += new_pos.x;

triangle.vertex1.y += new_pos.y;

triangle.vertex1.z += new_pos.z;

triangle.vertex2.x += new_pos.x;

triangle.vertex2.y += new_pos.y;

triangle.vertex2.z += new_pos.z;

triangle.vertex3.x += new_pos.x;

triangle.vertex3.y += new_pos.y;

triangle.vertex3.z += new_pos.z;

}

//Rodrig rotation

void rotate(const Point3d& _pos, Model_Triangle& triangle, float angle){

Vector3d z_axis(0.0, 0.0, 1.0);

Vector3d pos(_pos);

pos.normalize();

Vector3d rot_axis = z_axis * pos;

rot_axis.normalize();

float cos_theta = cos(angle* 3.14159265 / 180.0);

float sin_theta = sin(angle * 3.14159265 / 180.0);

Vector3d v1(triangle.vertex1);

Vector3d v2(triangle.vertex1);

Vector3d v3(triangle.vertex1);

Vector3d v4 = triangle.v3d_normal;

Vector3d kv_cross_v1 = rot_axis*v1;

Vector3d kv_cross_v2 = rot_axis*v2;

Vector3d kv_cross_v3 = rot_axis*v3;

Vector3d kv_cross_v4 = rot_axis*v4;

float kv_dot1 = rot_axis.dot(v1);

float kv_dot2 = rot_axis.dot(v2);

float kv_dot3 = rot_axis.dot(v3);

float kv_dot4 = rot_axis.dot(v4);

//Rotate with Rodriges formula

triangle.vertex1.x = v1.x * cos_theta + kv_cross_v1.x * sin_theta +

rot_axis.x * kv_dot1 * (1 - cos_theta);

triangle.vertex1.y = v1.y * cos_theta + kv_cross_v1.y * sin_theta +

rot_axis.y * kv_dot1 * (1 - cos_theta);

triangle.vertex1.z = v1.z * cos_theta + kv_cross_v1.z * sin_theta +

rot_axis.z * kv_dot1 * (1 - cos_theta);

triangle.vertex2.x = v2.x * cos_theta + kv_cross_v2.x * sin_theta +

rot_axis.x * kv_dot2 * (1 - cos_theta);

triangle.vertex2.y = v2.y * cos_theta + kv_cross_v2.y * sin_theta +

rot_axis.y * kv_dot2 * (1 - cos_theta);

triangle.vertex2.z = v2.z * cos_theta + kv_cross_v2.z * sin_theta +

rot_axis.z * kv_dot2 * (1 - cos_theta);

triangle.vertex3.x = v3.x * cos_theta + kv_cross_v3.x * sin_theta +

rot_axis.x * kv_dot3 * (1 - cos_theta);

triangle.vertex3.y = v3.y * cos_theta + kv_cross_v3.y * sin_theta +

rot_axis.y * kv_dot3 * (1 - cos_theta);

triangle.vertex3.z = v3.z * cos_theta + kv_cross_v3.z * sin_theta +

rot_axis.z * kv_dot3 * (1 - cos_theta);

triangle.v3d_normal.x = v4.x * cos_theta + kv_cross_v4.x * sin_theta +

rot_axis.x * kv_dot4 * (1 - cos_theta);

triangle.v3d_normal.y = v4.y * cos_theta + kv_cross_v4.y * sin_theta +

rot_axis.y * kv_dot4 * (1 - cos_theta);

triangle.v3d_normal.z = v4.z * cos_theta + kv_cross_v4.z * sin_theta +

rot_axis.z * kv_dot4 * (1 - cos_theta);

}

//LookAt transofrmation of scene in camera view

void cam_view(const Point3d& _camera_pos,

const Vector3d _camera_normal,

Model_Triangle& triangle) {

Vector3d camera_pos(_camera_pos);

Vector3d up(0.0, 1.0, 0.0);

if (abs(_camera_normal.y) == 1.0) {

up.x = 1.0;

up.y = 0.0;

up.z = 0.0;

}

//Vector3d camera_normal = camera_pos + _camera_normal;

//define cam coordinate system

Vector3d zaxis;

//zaxis = camera_pos - camera_normal;

zaxis.x = -_camera_normal.x;

zaxis.y = -_camera_normal.y;

zaxis.z = -_camera_normal.z;

Vector3d xaxis = up*zaxis;

Vector3d yaxis = zaxis*xaxis;

Model_Triangle tmp;

//Rotate, meaning that multiplicate by transposed

//and translate

tmp.vertex1.x = dot_product(xaxis, triangle.vertex1) - dot_product(xaxis,

camera_pos);

tmp.vertex1.y = dot_product(yaxis, triangle.vertex1) - dot_product(yaxis,

camera_pos);

tmp.vertex1.z = dot_product(zaxis, triangle.vertex1) - dot_product(zaxis,

camera_pos);

tmp.vertex2.x = dot_product(xaxis, triangle.vertex2) - dot_product(xaxis,

camera_pos);

tmp.vertex2.y = dot_product(yaxis, triangle.vertex2) - dot_product(yaxis,

camera_pos);

tmp.vertex2.z = dot_product(zaxis, triangle.vertex2) - dot_product(zaxis,

camera_pos);

tmp.vertex3.x = dot_product(xaxis, triangle.vertex3) - dot_product(xaxis,

camera_pos);

tmp.vertex3.y = dot_product(yaxis, triangle.vertex3) - dot_product(yaxis,

camera_pos);

tmp.vertex3.z = dot_product(zaxis, triangle.vertex3) - dot_product(zaxis,

camera_pos);

tmp.v3d_normal.x = dot_product(xaxis, triangle.v3d_normal) -

dot_product(xaxis, camera_pos);

tmp.v3d_normal.y = dot_product(yaxis, triangle.v3d_normal) -

dot_product(yaxis, camera_pos);

tmp.v3d_normal.z = dot_product(zaxis, triangle.v3d_normal) -

dot_product(zaxis, camera_pos);

triangle = tmp;

}

//ray-triangle intersection

bool isintersect(Model_Triangle& target, Vector3d& D, Vector3d &dist, int

n) {

float t = 0;

float u = 0;

float v = 0;

Vector3d t_vec;

float norm_t = 0;

Vector3d E1((target.vertex2.x - target.vertex1.x),

(target.vertex2.y - target.vertex1.y),

(target.vertex2.z - target.vertex1.z));

Vector3d E2((target.vertex3.x - target.vertex1.x),

(target.vertex3.y - target.vertex1.y),

(target.vertex3.z - target.vertex1.z));

Vector3d T((-target.vertex1.x),

(-target.vertex1.y),

(-target.vertex1.z));

Vector3d P = D * E2;

Vector3d Q = T * E1;

float Z = P.dot(E1);

t = Q.dot(E2) / Z;

u = P.dot(T) / Z;

v = Q.dot(D) / Z;

t_vec.x = t*D.x;

t_vec.y = t*D.y;

t_vec.z = t*D.z;

dist = t_vec;

bool boo;

if ((u >= 0 && v >= 0) && ((u + v) <= 1) ) {

//inlog << "true" << endl;

boo = true;

} else {

//inlog << "false" << endl;

boo = false; }

return boo;

}

//ray-triangle intersection

bool isintersectsun(Model_Triangle& target, Vector3d& D, vector<float>

&dist, int n) {

ofstream inlog;

inlog.open("intrsctn_log.txt", ios_base::out | ios_base::app);

float t = 0;

float u = 0;

float v = 0;

Vector3d t_vec;

float norm_t = 0;

Vector3d E1((target.vertex2.x - target.vertex1.x),

(target.vertex2.y - target.vertex1.y),

(target.vertex2.z - target.vertex1.z));

Vector3d E2((target.vertex3.x - target.vertex1.x),

(target.vertex3.y - target.vertex1.y),

(target.vertex3.z - target.vertex1.z));

Vector3d T((-target.vertex1.x),

(-target.vertex1.y),

(-target.vertex1.z));

Vector3d P = D * E2;

Vector3d Q = T * E1;

float Z = P.dot(E1);

t = Q.dot(E2) / Z;

u = P.dot(T) / Z;

v = Q.dot(D) / Z;

t_vec.x = t*D.x;

t_vec.y = t*D.y;

t_vec.y = t*D.y;

norm_t = t_vec.norm();

bool boo;

if ((u >= 0 && v >= 0) && ((u + v) < 1) && (t > 0)) {

boo = true;

}

else {

boo = false;

}

return boo;

}

Приложение 10

func.cpp

#define _CRT_SECURE_NO_WARNINGS

#include "FUNC.H"

#include <iostream>

#include <fstream>

#include <iomanip>

#include <vector>

#include <cmath>

#define ABCORR "NONE"

using namespace std;

const float PI = 1 / 3.1415;

const float sb_const = 5.670367e-08;

const float moon_T = 271.5;

const float sun_brightness = 2e+07; // Wt/m2/str

const float sun_area = 3e+18; //actually, half of the area in m2

float sun_flux = sun_area*sun_brightness;

//Axuillary

Vector3d Distance(Model_Triangle& _i, Model_Triangle& _j) {

std::vector<float> temp_dist;

temp_dist.push_back((_j.p3d_centr.x) - (_i.p3d_centr.x));

temp_dist.push_back((_j.p3d_centr.y) - (_i.p3d_centr.y));

temp_dist.push_back((_j.p3d_centr.z) - (_i.p3d_centr.z));

Vector3d v3d_dist(temp_dist);

return v3d_dist;

}

float Norm_v3d(Vector3d& _v3d) {

float temp = 0;

temp = pow(_v3d.x, 2.0)

+ pow(_v3d.y, 2.0)

+ pow(_v3d.z, 2.0);

return (sqrt(temp));

}

float dot_product(Vector3d& a, Vector3d& b) {

float dot = 0;

dot = (a.x)*(b.x) + (a.y)*(b.y) + (a.z)*(b.z);

return dot;

}

float dot_product(Vector3d& a, Point3d& b) {

float dot = 0;

dot = (a.x)*(b.x) + (a.y)*(b.y) + (a.z)*(b.z);

return dot;

}

int find_min(vector<float> t) {

float min = 1000000;

int min_num = 0;

for (int k = 0; k < t.size(); k++){

if (t[k] != 0){

if (t[k] < min) {

min = t[k];

min_num = k;

}

else { continue; }

}

}

return min_num;

}

float char_to_float1(string& char_buffer) {

float fl = 0;

char* temp = new char[char_buffer.length() + 1];

std::strcpy(temp, char_buffer.c_str());

fl = atof(temp);

delete[] temp;

return fl;

}

//main functions

//computing view factors of

float compute_SURF_alpha(

Vector3d& _i_normal,

Vector3d& _j_normal,

Vector3d _ij_distance,

float(*dot_product)(Vector3d&, Vector3d&),

float(*Norm_v3d)(Vector3d&),

float _j_Area) {

float tmp_alpha = 0.;

float tmp_dp_i = dot_product(_i_normal, _ij_distance);

float tmp_dp_j = dot_product(_j_normal, _ij_distance);

float tmp_dp_ij = dot_product(_j_normal, _i_normal);

float tmp_dist_norm = (Norm_v3d(_ij_distance));

//tmp_alpha = PI*tmp_dp_i*tmp_dp_j / (pow(tmp_dist_norm,

4.0))*_j_Area;

if ((tmp_dp_ij <= 0))

tmp_alpha = PI*abs((tmp_dp_i*tmp_dp_j / (pow(tmp_dist_norm,

4.0)))*_j_Area);

else tmp_alpha = 0;

return tmp_alpha;

}

//computing view factors to calculate direct sun flux

float compute_SURF_alpha_sun(

Vector3d& _i_normal,

Vector3d& _j_normal,

Vector3d _ij_distance,

float(*dot_product)(Vector3d&, Vector3d&),

float(*Norm_v3d)(Vector3d&),

float _j_Area) {

Vector3d _ij_distance_n(_ij_distance);

_ij_distance_n.normalize();

float tmp_alpha = 0.;

float tmp_dp_j = dot_product(_j_normal, _ij_distance_n);

float tmp_dp_ij = dot_product(_j_normal, _i_normal);

float tmp_dist_norm = (Norm_v3d(_ij_distance));

if ((tmp_dp_ij <= 0))

tmp_alpha = PI*abs((tmp_dp_j / (pow(tmp_dist_norm*1000,

2.0)))*_j_Area/10000

);

else tmp_alpha = 0;

return tmp_alpha;

}

//determ visabillity and make view factor matrix

void make_matrix(deque<Model_Triangle>& SCENE, vector<

vector<float> > &reflection_matrix) {

ofstream log;

log.open("ViewFactorMatrix.txt");

int i;

for (i = 0; i < SCENE.size(); i++) { //main camera loop

std::cout << i << " iteration starts" << endl;

// put cam in a centr of i triangle

Point3d Camera(SCENE[i].p3d_centr.x, SCENE[i].p3d_centr.y,

SCENE[i].p3d_centr.z);

//Zaxis

Vector3d Camera_normal = SCENE[i].v3d_normal;

std::cout << i << " camera setted" << endl;

//A que with projected triangles

deque<Model_Triangle> projected;

for (int j = 0; j < SCENE.size(); j++){

if (i != j) {

//Temporary triangle preform transformations to

Model_Triangle tmp = SCENE[j];

//Transform scene to Camera View

cam_view(Camera, Camera_normal, tmp);

projected.emplace_back(tmp);

}

else { continue; }

} //After this loop we have a que with projected triangles to a setted cam

position

std::cout << i << " transformation done" << endl;

//Iterate through triangles with intersection finding

for (int j = 0; j < projected.size(); j++){

//Store centr of a projected triangle of interest

Point3d centr(

(projected[j].vertex1.x + projected[j].vertex2.x + projected[j].vertex3.x)/3,

(projected[j].vertex1.y + projected[j].vertex2.y + projected[j].vertex3.y)/3,

(projected[j].vertex1.z + projected[j].vertex2.z + projected[j].vertex3.z)/3

);

//Normallized direction vector

Vector3d D(centr);

D.normalize();

//not normalized D vector

Vector3d Dnn(centr);

int scene_n;

//we want to know what SCENE triangle been hit

if (j < i) { scene_n = j; }

else if (j >= i) { scene_n = j + 1; }

//i is a dist, j is a num

bool flag = false;

float tn = 0;

if ((Dnn.z > 0) || (abs(Dnn.z) < 0.00001)) { continue; }

//because it means we are looking at the triangle from the back

for (int n = 0; n < projected.size(); n++) {

// now we are seeking for intersections and cmp alpha

if (n == j) { continue; }

//distance of a intersected

Vector3d t;

if (isintersect(projected[n], D, t, scene_n) == true) {

//std::cout << abs(t.z) << " " << abs(Dnn.z) << endl;

if ( abs(t.z) <= abs(Dnn.z) ) {

flag = true;

break;

}

}

} // projected triangles loop done

if (flag == false) {

//std::cout << "cmpiting alpha now" << endl;

reflection_matrix[i][scene_n] =

-0.073*

compute_SURF_alpha(

SCENE[i].v3d_normal,

SCENE[scene_n].v3d_normal,

Distance(SCENE[i], SCENE[scene_n]),

&dot_product, &Norm_v3d, SCENE[scene_n].Area()

);

}

} // j triangle loop done

/*

//close pixel loop

}

}

*/

std::cout << "Now next cam triangle ->" << endl;

} //Main camera loop done

//Print result

std::cout << "printing matrix..." << endl;

for (int ni = 0; ni < reflection_matrix.size(); ni++){

for (int nj = 0; nj < reflection_matrix[ni].size(); nj++){

log << setiosflags(ios::fixed | ios::showpoint) << setprecision(10) <<

setw(12) << reflection_matrix[ni][nj] << " ";

}

log << endl;

}

//end of function

}

//solve to determ reflection fluxes

void zeidel(deque<Model_Triangle>& SCENE, vector<vector<float> >&

Matrix,

vector<float>& x, int columns,

float(*char_to_float1)(string&),

bool(*converge)(vector<float>& xk, vector<float>& xkp, int n)) {

ofstream verify;

verify.open("verify.txt");

string char_buffer;

vector<float> alphas_row;

float float_buffer = 0;

std::cout << "columns: " << columns << endl;

//vector< vector < float > > Matrix;

vector<float> b;

float tmp_StefBolz = sb_const*pow(moon_T, 4.0); //Stephan-Boltzman

radiance law

for (int i = 0; i < columns; i++)

b.push_back(

SCENE[i].sun_direct_B +

SCENE[i].sun_reflected_B +

SCENE[i].moon_temperature_B +

SCENE[i].selfreflection_B

);

for (int i = 0; i < b.size(); i++)

verify << "b[" << i << "]:" << b[i] << endl;

vector<float> p(columns, 0);

vector<float>::iterator it;

do {

p = x;

for (int i = 0; i < Matrix.size(); i++) {

float var = 0;

for (int j = 0; j < i; j++) {

var += (Matrix[i][j] * x[j]);

}

for (int j = i + 1; j < Matrix.size(); j++){

var += (Matrix[i][j] * p[j]);

}

x[i] = (b[i] - var);

}

}

while (!converge(x, p, Matrix.size()));

}

//convergence of zeidel

bool converge(vector<float>& xk, vector<float>& xkp, int n)

{

float norm = 0;

for (int i = 0; i < n; i++) { norm += (xk[i] - xkp[i])*(xk[i] - xkp[i]); }

std::cout << "Norm: " << norm << endl;

if (sqrt(norm) >(0.0001)) { return false; }

else { return true; }

}

//Compute fluxes

void cmp_fluxes_on_surf(deque<Model_Triangle>& SURF,

SpiceDouble et ) { //we are passing time of interest

std::ofstream out;

out.open("out213.txt", std::ios_base::out | std::ios_base::app | std::ios::right);

//compute flux from moon temperature

for (int i = 0; i < SURF.size(); i++) {

SURF[i].moon_temperature_B = sb_const*pow(moon_T, 4.0); // Wt/m2

}

//Determ occultation state to compute fluxes from sun

if(is_oculted_moon(et) == false) {

//Inputs are 80 ю.ш. 48 з.д

SpiceDouble latitude = (-80); //should be in radians when passing to a func

SpiceDouble longtitude = (-48); //should be in radians when passing to a

func

//output rec coord

SpiceDouble landing_pt[3];

//radiuses of 3axial elipsoid of the moon

SpiceDouble radiuses[3];

SpiceInt dim = 3;

//SpiceInt i;

//outpusts of ilum routine

SpiceDouble emissn;

SpiceDouble srfvec[3];

SpiceDouble phase;

SpiceDouble solar; //we need this. it's an angle from surf normal to sun

position vector

SpiceDouble trgepc;

//time strings

SpiceChar step_time[51];

SpiceChar time[51];

//Convert to radians

latitude *= rpd_c();

longtitude *= rpd_c();

//get moon radii

bodvrd_c("MOON", "RADII", 3, &dim, radiuses);

//convert to cartesian

latrec_c(radiuses[0], longtitude, latitude, landing_pt);

//make vector3d from landing_pt

Point3d land_pt(landing_pt);

//the angle from zaxis to latrec vector

// its 170 degrees because rotation is anticlockwise

float lat_grad = -(dpr_c()*latitude - 90);

//here we compute solar angle. other stuff is not important

ilumin_c("Ellipsoid", "MOON", et, "IAU_MOON", "CN+S", "SUN",

landing_pt,

&trgepc, srfvec, &phase, &solar, &emissn);

et2utc_c(et, "C", 0, 22, step_time);

out << step_time << " solar: " << dpr_c()*solar << endl;

//cout << "surf solar: " << dpr_c()*solar << endl;

SpiceDouble distance_to_p = vnorm_c(srfvec)*1000; //distance in

killometers converted to meters

Vector3d surf_vec(srfvec);

Vector3d surf_vec_n(srfvec);

surf_vec_n.normalize();

deque<Model_Triangle> tmp_SURF;

for (int i = 0; i < SURF.size(); i++){

tmp_SURF.emplace_back(SURF[i]);

}

for (int i = 0; i < tmp_SURF.size(); i++){

rotate(land_pt, tmp_SURF[i], lat_grad);

}

//compute B for all surface triangles

for (int i = 0; i < tmp_SURF.size(); i++) {

float tmp_alpha = 0.;

Vector3d temp(tmp_SURF[i].v3d_normal.x, tmp_SURF[i].v3d_normal.z,

tmp_SURF[i].v3d_normal.y);

if (temp.dot(surf_vec_n) <= 0){

tmp_alpha =

abs(

(PI*cos(solar) / (pow(distance_to_p, 2.0)))

*(SURF[i].Area()/10000) //area in sm converted to meters

);

}

else { tmp_alpha = 0; }

//set radiosity on each triangle in Wt/m2

SURF[i].sun_direct_B = tmp_alpha*(1 - SURF[i].p)*sun_flux;

SURF[i].sun_reflected_B = tmp_alpha*(SURF[i].p)*sun_flux;

out << "tmp_alpha: " << tmp_alpha << endl

<< "sun_direct_B: " << SURF[i].sun_direct_B << endl

<< "sun_reflected_B: " << SURF[i].sun_reflected_B << endl

<< "moon_temperature_B: " << SURF[i].moon_temperature_B << endl;

}

} else { cout << "Occulted" << endl; }

out.close();

}

void cmp_direct_flux_on_box(deque<Model_Triangle>& SCENE,

SpiceDouble et){

if (is_oculted_moon(et) == false) {

//Inputs are 80 ю.ш. 48 з.д

SpiceDouble latitude = (-80);

SpiceDouble longtitude = (-48);

//output rec coord

SpiceDouble landing_pt[3];

//Moon triaxial elipsoid radiuses

SpiceDouble radiuses[3];

SpiceInt dim = 3;

//SpiceInt i;

SpiceDouble lt;

//angle and sun-pt vec

SpiceDouble solar;

SpiceDouble sun_pt_vec[6];

//time strings

SpiceChar step_time[51];

SpiceChar time[51];

// Copy triangles to preform action to.

// We change Y and Z coordinates because in model space Z it's a depth

// and we change units from mm to km

deque<Model_Triangle> BOX;

for (int i = 0; i < SCENE.size(); i++) {

BOX.push_back(SCENE[i]);

float tmp_y = BOX[i].vertex1.y;

float tmp_z = BOX[i].vertex1.z;

BOX[i].vertex1.y = tmp_z / 100000;

BOX[i].vertex1.z = tmp_y / 100000;

tmp_y = BOX[i].vertex2.y;

tmp_z = BOX[i].vertex2.z;

BOX[i].vertex2.y = tmp_z / 100000;

BOX[i].vertex2.z = tmp_y / 100000;

tmp_y = BOX[i].vertex3.y;

tmp_z = BOX[i].vertex3.z;

BOX[i].vertex3.y = tmp_z / 100000;

BOX[i].vertex3.z = tmp_y / 100000;

tmp_y = BOX[i].v3d_normal.y;

tmp_z = BOX[i].v3d_normal.z;

BOX[i].v3d_normal.y = tmp_z;

BOX[i].v3d_normal.z = tmp_y;

}

//Convert to radians

latitude *= rpd_c();

longtitude *= rpd_c();

//get moon radii

bodvrd_c("MOON", "RADII", 3, &dim, radiuses);

//convert to cartesian

latrec_c(radiuses[0], longtitude, latitude, landing_pt);

//make vector3d from landing_pt

Point3d land_pt(landing_pt);

//the angle from zaxis to latrec vector

// its 170 degrees because rotation is anticlockwise

float lat_grad = -(dpr_c()*latitude - 90);

//now compute angles

std::ofstream out;

out.open("out.txt", std::ios_base::out | std::ios_base::app | std::ios::right);

//get vector sun-lnd_pt

spkcpt_c(landing_pt, "MOON", "IAU_MOON", et, "IAU_MOON", "OBSERVER", ABCORR,

"SUN", sun_pt_vec, &lt);

et2utc_c(et, "C", 0, 22, step_time);

float sun_pt_3[3];

sun_pt_3[0] = sun_pt_vec[0];

sun_pt_3[1] = sun_pt_vec[1];

sun_pt_3[2] = sun_pt_vec[2];

//length of sun-point vector

SpiceDouble distance_to_p = vnorm_c(sun_pt_vec) * 1000; //in meters

Vector3d sun_pt_v3d(sun_pt_3);

Vector3d sun_pt_v3d_n(sun_pt_3);

sun_pt_v3d_n.normalize();

Vector3d landing_pt_v3d(landing_pt);

//compute B for all BOX triangles

for (int i = 0; i < BOX.size(); i++) {

//transform to IAU_MOON frame

rotate(land_pt, BOX[i], lat_grad);

//translate(land_pt, BOX[i]);

float tmp_alpha = 0.;

//compute angle between normal and sun_pt_v3d in radians

solar = acos(BOX[i].v3d_normal.dot(sun_pt_v3d_n));

//convert to radians

if (BOX[i].v3d_normal.dot(sun_pt_v3d_n) <= 0) {

tmp_alpha = abs((PI*cos(solar) / (pow(distance_to_p,

2.0)))//*(SCENE[i].Area() / 10000))

);

//distance units is m, area converted from sm to m

}

else tmp_alpha = 0;

//Put computed values in BOX SCENE triangles instead of copyed

SCENE[i].sun_direct_B = tmp_alpha*(1 - SCENE[i].p)*sun_flux; // radiosity in Wt/m2

}

} else { cout << "Occulted" << endl; }

}

void cmp_model_flux_from_surf(deque<Model_Triangle>& SCENE,

deque<Model_Triangle>& SURF) {

for (int i = 0; i < SCENE.size(); i++) {

float alpha = 0;

for (int j = 0; j < SURF.size(); j++){

alpha = compute_SURF_alpha(SURF[j].v3d_normal,

SCENE[i].v3d_normal,

Distance(SURF[j], SCENE[i]),

&dot_product, &Norm_v3d, SCENE[i].Area());

SCENE[i].sun_reflected_B += alpha*SURF[j].sun_reflected_B;//

*SURF[i].Area() / 100;

SCENE[i].moon_temperature_B += alpha*SURF[j].moon_temperature_B;//

*SURF[i].Area() / 100;

}

}

}

void cmp_direct_flux_on_model(deque<Model_Triangle>& SCENE_,

SpiceDouble et) {

//First determ occultation of the landing point

if (is_oculted_moon(et) == false) {

//copy to tmp SCENE to rotate it

deque <Model_Triangle> SCENE;

for (int i = 0; i < SCENE_.size(); i++)

SCENE.emplace_back(SCENE_[i]);

//Inputs are 80 ю.ш. 48 з.д

SpiceDouble latitude = (-80);

SpiceDouble longtitude = (-48);

//output rec coord

SpiceDouble landing_pt[3];

SpiceDouble radiuses[3];

SpiceInt dim = 3;

//SpiceInt i;

SpiceDouble lt;

SpiceDouble emissn;

SpiceDouble srfvec[3];

SpiceDouble phase;

SpiceDouble solar;

SpiceDouble trgepc;

SpiceDouble sun_pt_vec[6];

SpiceChar step_time[51];

SpiceChar time[51];

//str2et_c(time2, &et2);

//str2et_c(time, &et);

//Convert to radians

latitude *= rpd_c();

longtitude *= rpd_c();

//get moon radii

bodvrd_c("MOON", "RADII", 3, &dim, radiuses);

//convert to cartesian

latrec_c(radiuses[0], longtitude, latitude, landing_pt);

//make vector3d from landing_pt

Point3d land_pt(landing_pt);

float lat_grad = -(dpr_c()*latitude - 90);

//std::cout << "latitude: " << lat_grad << endl;

//Set model to IAU_MOON Frame on the landing point

for (int i = 0; i < SCENE.size(); i++) {

rotate(land_pt, SCENE[i], lat_grad);

translate(land_pt, SCENE[i]);

}

// now the scene is setted to compute angles

std::ofstream out;

out.open("out.txt", std::ios_base::out | std::ios_base::app | std::ios::right);

//get vector lnd_pt - sun

spkcpt_c(landing_pt, "MOON", "IAU_MOON", et, "IAU_MOON",

"OBSERVER", ABCORR,

"SUN", sun_pt_vec, &lt);

SpiceDouble sun_moon[3];

//get sun position in IAU_MOON

spkpos_c("SUN", et, "IAU_MOON", ABCORR, "MOON", sun_moon, &lt);

et2utc_c(et, "C", 0, 22, step_time);

out << step_time << " sun_pt_vec: " << endl

<< sun_pt_vec[0] << endl

<< sun_pt_vec[1] << endl

<< sun_pt_vec[2] << endl;

float sun_pt_3[3];

sun_pt_3[0] = sun_pt_vec[0];

sun_pt_3[1] = sun_pt_vec[1];

sun_pt_3[2] = sun_pt_vec[2];

SpiceDouble distance_to_p = vnorm_c(sun_pt_3);

float sun_moon_distance = vnorm_c(sun_moon);

Vector3d sun_pt_v3d(sun_pt_3);

Vector3d sun_moon_v3d(sun_moon); //in km

//Now we need to determ intersected triangles

// put cam in SUN

Point3d Camera(sun_moon);

//Zaxis

Vector3d Camera_normal(sun_pt_3);

Camera_normal.normalize();

//A que with projected triangles

deque<Model_Triangle> projected;

for (int j = 0; j < SCENE.size(); j++){

//Temporary triangle preform transformations to

Model_Triangle tmp = SCENE[j];

//Transform scene to Camera View

cam_view(Camera,

Camera_normal,

tmp);

projected.emplace_back(tmp);

}

//After this loop we have a que with projected triangles to a setted cam

position

std::cout << "sun lookAt transformation done" << endl;

for (int j = 0; j < projected.size(); j++){

//Store centr of a projected triangle of interest

Point3d centr(

(projected[j].vertex1.x + projected[j].vertex2.x + projected[j].vertex3.x) / 3,

(projected[j].vertex1.y + projected[j].vertex2.y + projected[j].vertex3.y) / 3,

(projected[j].vertex1.z + projected[j].vertex2.z + projected[j].vertex3.z) / 3

);

//Normallized direction vector

Vector3d D(centr);

D.normalize();

//not normalized D vector

Vector3d Dnn(centr);

bool flag = false;

if ((Dnn.z > 0) || (abs(Dnn.z) < 0.00001)) { continue; } //because it means

we are looking at the triangle from the back

for (int n = 0; n < projected.size(); n++) { // now we are seeking for

intersections and cmp alpha

if (n == j) { continue; }

//distance of a intersected

Vector3d t;

if (isintersect(projected[n], D, t, n) == true) {

if (abs(t.z) <= abs(Dnn.z)) {

flag = true;

break;

}

}

} // projected triangles loop done

if (flag == false) {

//std::cout << "cmpiting alpha now" << endl;

SCENE_[j].sun_direct_B =

0.93*

compute_SURF_alpha_sun(

Dnn,

SCENE[j].v3d_normal,

Dnn,

&dot_product, &Norm_v3d, SCENE_[j].Area())*sun_flux;

}

}// j triangle loop done

}

else { std::cout << "Landing point is in the moon shadow" << endl; }

//end of function

}

bool is_oculted_moon(SpiceDouble et) {

SpiceDouble avg;

SpiceDouble r;

SpiceDouble s[2][6];

SpiceDouble sep;

SpiceBoolean found;

SpiceInt j;

SpiceInt k;

SpiceInt t[2];

//Inputs are 80 ю.ш. 48 з.д

SpiceDouble latitude = (-80);

SpiceDouble longtitude = (-48);

//output rec coord

SpiceDouble landing_pt[3];

SpiceDouble radiuses[3];

SpiceInt dim = 3;

SpiceDouble lt;

SpiceDouble emissn;

SpiceDouble srfvec[3];

SpiceDouble phase;

SpiceDouble solar;

SpiceDouble trgepc;

SpiceDouble sun_pt_vec[6];

SpiceChar step_time[51];

SpiceChar time[51];

//Convert to radians

latitude *= rpd_c();

longtitude *= rpd_c();

//get moon radii

bodvrd_c("MOON", "RADII", 3, &dim, radiuses);

//convert to cartesian

latrec_c(radiuses[0], longtitude, latitude, landing_pt);

//get vector lnd_pt - sun

spkcpt_c(landing_pt, "MOON", "IAU_MOON", et, "IAU_MOON",

"OBSERVER", ABCORR,

"SUN", sun_pt_vec, &lt);

//make vector3d from landing_pt

Point3d land_pt(landing_pt);

Vector3d ldpt(landing_pt);

ldpt.norm();

//check occultation

//Get the ID codes associated with the targets

bodn2c_c("earth", &t[0], &found);

bodn2c_c("sun", &t[1], &found);

avg = sumad_c(radiuses, 3) / 3.0;

r = asin(avg / ldpt.norm());

/*

Determine the separation between the two bodies. If the

separation between the centers is greater than the sum of

the apparent radii, then the target bodies are clear of

each other.

Function vsep_c returns the angle of separation between

two three-vectors.

*/

sep = vsep_c(sun_pt_vec, landing_pt) - r;

string strng;

if (sep > 0.)

{

cout << "Clear" << endl;

return false;

}

else

{

return true;

}

}

void find_sun() {

SpiceDouble avg;

SpiceDouble r;

SpiceDouble s[2][6];

SpiceDouble sep;

SpiceBoolean found;

SpiceInt j;

SpiceInt k;

SpiceInt t[2];

//Inputs are 80 ю.ш. 48 з.д

SpiceDouble latitude = (-80);

SpiceDouble longtitude = (-48);

//output rec coord

SpiceDouble landing_pt[3];

SpiceDouble radiuses[3];

SpiceInt dim = 3;

SpiceDouble lt;

SpiceInt i;

SpiceDouble emissn;

SpiceDouble srfvec[3];

SpiceDouble phase;

SpiceDouble solar;

SpiceDouble trgepc;

SpiceDouble sun_pt_vec[6];

SpiceChar step_time[51];

SpiceChar time[51];

SpiceChar time1[51];

SpiceChar time2[51];

SpiceDouble et;

SpiceDouble et1;

SpiceDouble et2;

#define ET0 0.0

/*

Use a time step of 1 hour; look up 100 states.

*/

#define STEP 3600.0

#define MAXITR 1000

#define TIMLEN 51

//Specify epoch

prompt_c("Epoch from? ", 51, time1);

prompt_c("Epoch to? ", 51, time2);

str2et_c(time1, &et1);

str2et_c(time2, &et2);

; i < MAXITR + 1; i++)

{

et = et1 + (i*(et2 - et1) / MAXITR);

//get moon radii

bodvrd_c("MOON", "RADII", 3, &dim, radiuses);

//convert to cartesian

latrec_c(radiuses[0], longtitude, latitude, landing_pt);

//get vector lnd_pt - sun

spkcpt_c(landing_pt, "MOON", "IAU_MOON", et, "IAU_MOON",

"OBSERVER", ABCORR,

"SUN", sun_pt_vec, &lt);

//make vector3d from landing_pt

Point3d land_pt(landing_pt);

Vector3d ldpt(landing_pt);

ldpt.norm();

//check occultation

//Get the ID codes associated with the targets

bodn2c_c("earth", &t[0], &found);

bodn2c_c("sun", &t[1], &found);

avg = sumad_c(radiuses, 3) / 3.0;

r = asin(avg / ldpt.norm());

/*

Determine the separation between the two bodies. If the

separation between the centers is greater than the sum of

the apparent radii, then the target bodies are clear of

each other.

Function vsep_c returns the angle of separation between

two three-vectors.

*/

sep = vsep_c(sun_pt_vec, landing_pt) - r;

string strng;

et2utc_c(et, "C", 0, 22, step_time);

if (sep > 0.)

{

cout << step_time << " Clear" << endl;

Размещено на Allbest.ru