Расчет посадки космического аппарата "Луна-Глоб" на Луну

//Нахождение длины вектора по координатам.

double length(double * vect)

{double l;

l = sqrt(pow(vect[0], 2) + pow(vect[1], 2) + pow(vect[2], 2));

return l;}

//Нахождение угла между векторами

double angle(double *v1, double *v2)

{double ang;

ang = acos((v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]) / (length(v1)*length(v2)));

return ang;}

//Расчет угла места КА над горизонтом относительно точки point

//point[3] - массив координат x y z точки на поверхности Луны

//KA[3] - массив координат x y z КА

double elevation_angle(double * point, double * KA)

{double Rl,Ra,r, betta, b, gamma, g, alpha;

double R[3];

double rtd = 57.295779513082320876798161804285;

Rl = length(point);//длина вектора Rl

Ra = length(KA);//длина вектора Ra

//Вектор R

R[0] = KA[0] - point[0];

R[1] = KA[1] - point[1];

R[2] = KA[2] - point[2];

r = length(R);//длина вектора R

betta = angle(point, KA);//угол между векторами Rl и Ra

gamma = angle(R, KA);//угол между векторами R и Ra

//Перевод из радиан в градусы.

b = betta*rtd;

g = gamma*rtd;

alpha = 90 - g - b;

return alpha;}

//умножение строки на матрицу

double * str_on_matrix(double m1[3], double m2[3][3])

{double * result = (double *)malloc(3 * sizeof(double));

result[0] = m1[0] * m2[0][0] + m1[1] * m2[1][0] + m1[2] * m2[2][0];

result[1] = m1[0] * m2[0][1] + m1[1] * m2[1][1] + m1[2] * m2[2][1];

result[2] = m1[0] * m2[0][2] + m1[1] * m2[1][2] + m1[2] * m2[2][2];

return result;}

// Перевод из СК SG в СК J2000.

double * SG_to_J(double Rsg[], double dt)

{double * Rj = (double *)malloc(3 * sizeof(double));

int i, j;//i-строка, j-столбец

double T, alph, bet, w, e[13], M1[3][3], M2[3][3];

int p = 10;

double dtr = 0.017453292519943295769236905555556; // коэффициент перевода градусов в радианы

double gaml = 13.17635815; // угл. скорость вращения Луны в градусах.

dt = dt - 2415020.0;

e[0] = (125.045 - 0.0529921*dt)*dtr;

e[1] = (250.089 - 0.1059842*dt)*dtr;

e[2] = (260.008 + 13.0120009*dt)*dtr;

e[3] = (176.625 + 13.3407154*dt)*dtr;

e[4] = (357.529 + 0.9856003*dt)*dtr;

e[5] = (311.589 + 26.4057084*dt)*dtr;

e[6] = (134.963 + 13.0649930*dt)*dtr;

e[7] = (276.617 + 0.3287146*dt)*dtr;

e[8] = (34.226 + 1.7484877*dt)*dtr;

e[9] = (15.134 - 0.1589763*dt)*dtr;

e[10] = (119.743 + 0.0036096*dt)*dtr;

e[11] = (239.961 + 0.1643573*dt)*dtr;

e[12] = (25.053 + 12.9590088*dt)*dtr;

T = dt / 36525;

alph = (269.9949 + 0.0031*T - 3.8787*sin(e[0]) - 0.1204*sin(e[1]) + 0.07*sin(e[2]) - 0.0172*sin(e[3]) + 0.0072*sin(e[5]) - 0.0052*sin(e[9]) - 0.0043*sin(e[12]))*dtr;

bet = (66.5392 + 0.013*T + 1.5419*cos(e[0]) + 0.0239*cos(e[1]) - 0.0278*cos(e[2]) + 0.0068*cos(e[3]) - 0.0029*cos(e[5]) + 0.0009*cos(e[6]) + 0.0008*cos(e[9]) - 0.0009*cos(e[12]))*dtr;

w = (38.3213 + gaml*dt + gaml + dt - 1.4*(p, (-12))*pow(dt, 2) + 3.561*sin(e[0]) + 0.1208*sin(e[1]) - 0.0642*sin(e[2]) + 0.0158*sin(e[3]) + 0.0252*sin(e[4]) - 0.0066*sin(e[5]) - 0.0047*sin(e[6]) - 0.0046*sin(e[7]) + 0.0028*sin(e[8]) + 0.0052*sin(e[9]) + 0.004*sin(e[10]) + 0.0019*sin(e[11]) - 0.0044*sin(e[12]))*dtr;

M1[0][0] = -sin(alph)*cos(w) - cos(alph)*sin(bet)*sin(w);

M1[0][1] = cos(alph)*cos(w) - sin(alph)*sin(bet)*sin(w);

M1[0][2] = sin(w)*cos(bet);

M1[1][0] = sin(alph)*sin(w) - cos(alph)*sin(bet)*cos(w);

M1[1][1] = -cos(alph)*sin(w) - sin(alph)*sin(bet)*cos(w);

M1[1][2] = cos(w)*cos(bet);

M1[2][0] = cos(alph)*cos(bet);

M1[2][1] = sin(alph)*cos(bet);

M1[2][2] = sin(bet);

//Траспонирование матрицы

for (i = 0; i < 3; i++)

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

{M2[i][j] = M1[j][i];}

}

Rj = str_on_matrix(Rsg, M2);

return Rj;}

int read_line(FILE * STREAM, char *line){

int i;

if (fgets(line, 256, STREAM) == NULL)

return (-1);

i = strlen(line);

while (i>0 && (line[i - 1] == '\r' || line[i - 1] == '\n')){

i--; line[i] = '\0';}

return i;}

/*Перевод из сферических в декартовы*/

void S_to_D(double sph[], double dec[])

{//широта долгота радиус

dec[0] = sph[2] * sin(sph[0])*sin(sph[1]);

dec[1] = sph[2] * cos(sph[0])*sin(sph[1]);

dec[2] = sph[2] * cos(sph[1]);}

int main()

{

setlocale(LC_ALL, "RUSSIAN");

int sch = 0, Nstr = 0;// sch-счетчик для case; Nstr-номер строки

char string[300], *pch;

//Объявление переменных

double R_Moon = 1738.57;// радиус Луны

double point[3] = { 175, 74, R_Moon };// точка на поверхности Луны

double decpoint[3] = { 0 };

double *p;

double x_sc, y_sc, z_sc;//координаты КА в СК J200

int year, month, day, hour, min, sec;

double KA[3];

char line[1024];

double JD = 0;//Юлианская дата

double ugol = 0;//угол места

double moon[3];

int npl = 10;

double Rl, Ra;//радиус Луны и расстояние до КА

double * shir, * dolg;//диапозоны по широте и долготе

char flag[10];//флаг. Если КА виден, пишет "да", если нет - "нет"

FILE * f;//f-файл, куда будем записывать

if ((f = fopen("D:\\image\\Катя\\Documents\\Visual Studio 2013\\Projects\\vidimost_bes_reliefa\\угол_места.txt", "wt")) == 0)

{printf("Ошибка!");

return 0;}

fprintf(f, " \t\t\t\t point \t\t\t\t\t\t\t\t\t\t KA J200 \t\t\t\t\t\t time \t\t\t\t угол места \t\t виден? \n");

fprintf(f, " point x \t point y \t\t point z \t\t KA J2000 x \t KA J2000 y \t KA J200 z \n");

//Переводим из сферических в декартовы

S_to_D(point, decpoint);

FILE * o;//o-файл, откуда читаемл

if ((o = fopen("D:\\image\\Катя\\Documents\\Visual Studio 2013\\Projects\\Проект 8(с файлами)\\ОРБИТА.txt", "rt")) == 0)

{printf("Ошибка!");

return 0;}

fgets(string, 300, o);

while (!feof(o))

{fgets(string, 300, o);

if (string != NULL)

{pch = strtok(string, " ");

KA[0] = atof(pch);

while (pch != NULL){

if (sch<5) pch = strtok(NULL, " ");

if ((sch>5) && (sch<9)) pch = strtok(NULL, " .");

if (sch >= 9) pch = strtok(NULL, ":");

switch (sch) {

case 0: KA[1] = atof(pch); break;

case 1: KA[2] = atof(pch); break;

case 6: day = atof(pch); break;

case 7: month = atof(pch); break;

case 8: year = atof(pch); break;

case 9: hour = atof(pch); break;

case 10: min = atof(pch); break;

case 11: sec = atof(pch); break;}

sch++;}

sch = 0;

//Переводим из UTC в Юлианскую дату

JD = DateTimeToJulDate(year, month, day, hour, min, sec);

//Переводим из SG в СК J2000

p = SG_to_J(decpoint, JD);

//Находим угол места

ugol = elevation_angle(p, KA);

//Определяем виден КА или нет.

if (ugol > 0)

{flag[0] = 'д';

flag[1] = 'а';

flag[2] = ' ';}

else

{flag[0] = 'н';

flag[1] = 'е';

flag[2] = 'т';}

fprintf(f, "%lf\t %lf\t %lf\t %15.5lf %15.5lf %15.5lf\t\t %d.%d.%d %d:%d:%d\t %15.5lf \t %5.5c%c%c \n", point[0], point[1], point[2], KA[0], KA[1], KA[2], year, month, day, hour, min, sec, ugol, flag[0], flag[1], flag[2]);

Nstr++;}}

fclose(f);

fclose(o);

return 0;}

Файл «main.cpp».

/*нахождение угла места КА и угла места высот и максимального времени виденья для заданной точки*/

#define _CRT_SECURE_NO_DEPRECATE

#include <time.h>

#include <stdio.h>

#include <cstdio>

#include <conio.h>

#include <locale.h>

#include <math.h>

#include <stdlib.h>

#include <string.h>

#include "orbita.h"

#include "moon.h"

#include "dop.h"

#define FILE_ORBIT "D:\\image\\Катя\\Documents\\Visual Studio 2013\\Projects\\vidimost_s_reliefom\\ОРБИТА.txt"

#define FILE_MOON "D:\\image\\Катя\\Desktop\\ДИПЛОМ\\Moon\\moon.txt"

#define FILE_ANGLES "D:\\image\\Катя\\Documents\\Visual Studio 2013\\Projects\\vidimost_s_reliefom\\углы.csv"

//"D:\\image\\Катя\\Documents\\Visual Studio 2013\\Projects\\vidimost_s_reliefom\\видимость.txt"

#define MoonRadius 1738.57

#define AngleStep 0.5e0

int main()

{int i, j, k;

orbita_point*KA;

int NKA;

int iO; // index of Orbit point

moon_height*MoonHeight;

double*angle_KA;

double*angle_Moon;

//Объявление переменных

double point_S[3] = { 65, 87, MoonRadius };// точка на поверхности Луны P'(lat,lon,R_Moon)

double point_J[3];

double MoonCenterSG[3];

double MoonCenterJ[3];

double vecOP[3],vecOK[3],vecN[3],vecOS[3];

double pointS_J[3],pointS_S[3];

double MatrixRotation[3][3];

double AngleNStep;

double a;

orbit_read(FILE_ORBIT, &KA, &NKA);

moon_height_read(FILE_MOON, 60, 90, 0, 360, -15520, 15520, &MoonHeight);

MoonCenterSG[0] = 0.e0;

MoonCenterSG[1] = 0.e0;

MoonCenterSG[2] = 0.e0;

angle_KA = (double*)calloc(NKA, sizeof(double));

angle_Moon = (double*)calloc(NKA, sizeof(double));

for (iO = 0; iO < NKA; iO++){

point_S[2] = MoonRadius + moon_height_get(MoonHeight, point_S[0], point_S[1]);

S_to_J(point_S, point_J, KA[iO].time);

SG_to_J(MoonCenterSG, MoonCenterJ, KA[iO].time);

for (k = 0; k < 3; k++) vecOP[k] = point_J[k] - MoonCenterJ[k];

for (k = 0; k < 3; k++) vecOK[k] = KA[iO].coords[k] - MoonCenterJ[k];

if (fabs(angle(vecOP, vecOK)) < 1.e-10){

//TODO видно, угол 90

angle_KA[iO] = 90.e0;

angle_Moon[iO] = 0.e0;

continue;}

angle_KA[iO] = elevation_angle(vecOP, vecOK);

mvector(vecOP, vecOK, vecN);

vec_norm(vecN);

AngleNStep = (int)(angle(vecOP, vecOK) / AngleStep + 1);

for (i = 1; i <= AngleNStep; i++){

matrix_rotation(vecN, i*AngleStep, MatrixRotation);

vec_on_matrix(vecOP, MatrixRotation, vecOS);

for (k = 0; k < 3; k++) pointS_J[k] = MoonCenterJ[k] + vecOS[k];

J_to_S(pointS_J, pointS_S,KA[iO].time);

pointS_S[2] = MoonRadius + moon_height_get(MoonHeight, pointS_S[0], pointS_S[1]);

S_to_J(pointS_S, pointS_J, KA[iO].time);

for (k = 0; k < 3; k++) vecOS[k] = pointS_J[k] - MoonCenterJ[k];

a = elevation_angle(vecOP, vecOS);

if (i == 1){

angle_Moon[iO] = a;

}else{

if (angle_Moon[iO] < a) angle_Moon[iO] = a;}}

printf("i=%d\tangle_KA=%g\tangle_Moon=%g\n", iO, angle_KA[iO], angle_Moon[iO]);}

setlocale(LC_ALL, "RUSSIAN");

FILE * fd;

if ((fd = fopen(FILE_ANGLES, "wt")) == 0){

printf("Can't open file Angels to write '%s'\n",FILE_ANGLES);

getch();

return 0;}

fprintf(fd, "Итерация;Дата;Время;Угол КА;Угол рельефа\n");

for (i = 0; i < NKA; i++){

fprintf(fd, "%d;%d-%d-%d;%d:%d:%d;%g;%g\n", i,

KA[i].time[0], KA[i].time[1], KA[i].time[2],

KA[i].time[3], KA[i].time[4], KA[i].time[5],

angle_KA[i], angle_Moon[i]);}

fclose(fd);

double TimeMax,Time,TimeStart,TimeFinish;

TimeMax = 0.e0;;

TimeStart = -1.e0;

for (i = 0; i < NKA; i++){

Time = DateTimeToJulDate(KA[i].time[0], KA[i].time[1], KA[i].time[2], KA[i].time[3], KA[i].time[4], KA[i].time[5]);

if (angle_KA[i] > angle_Moon[i]){

if (TimeStart < 0){

TimeStart = Time;

TimeFinish = Time;

}else{

TimeFinish = Time;}

}else{

if (TimeStart>0){

if (TimeFinish - TimeStart > TimeMax)

TimeMax = TimeFinish - TimeStart;

}else{//do nothing}

}}

if (TimeStart > 0){

if (TimeFinish - TimeStart > TimeMax)

TimeMax = TimeFinish - TimeStart;}

printf("Time max = %g days = %g sec\n", TimeMax,TimeMax*24*60*60);

getch();

return 0;}

Файл с вспомогательными программами «dop.cpp».

#define _CRT_SECURE_NO_DEPRECATE

#include <time.h>

#include <stdio.h>

#include <cstdio>

#include <conio.h>

#include <locale.h>

#include <math.h>

#include <stdlib.h>

#include <string.h>

#include "dop.h"

typedef struct{

double coord[3];

double JDtime;

}vector;

typedef struct{

double teta, fi, r;

//Получение Юлианской даты по Дате и Времени

double DateTimeToJulDate(int Year, int Month, int Day, int Hour, int Min, double Sec){ long MjdMidnight;

int b;

if (Month <= 2){ Month += 12; --Year; }

if ((10000L * Year + 100L * Month + Day) <= 15821004L)

b = -2 + ((Year + 4716) / 4) - 1179;

else

b = (Year / 400) - (Year / 100) + (Year / 4);

MjdMidnight = 365L * Year - 679004 + b + int(30.6001*(Month + 1)) + Day;

return MjdMidnight + ((Hour + Min / 60.0 + Sec / 3600.0) / 24.0) + 2400000.5;}

//Нахождение длины вектора по координатам.

double length(double * vect)

{double l;

l = sqrt(pow(vect[0], 2) + pow(vect[1], 2) + pow(vect[2], 2));

return l;}

//Скалярное произведение векторов

double skalyr(double *v1, double *v2){

return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];}

//Нахождение угла между векторами

double angle(double *v1, double *v2){

return RTD*acos(skalyr(v1, v2) / (length(v1)*length(v2)));}

//Векторное произведение векторов

void mvector(double v1[3], double v2[3],double v3[3]){

v3[0] = v1[1] * v2[2] - v1[2] * v2[1];

v3[1] = v1[2] * v2[0] - v1[0] * v2[2];

v3[2] = v1[0] * v2[1] - v1[1] * v2[0];}

//Нормировать вектор

void vec_norm(double v[3]){

double len;

len = length(v);

v[0] = v[0] / len;

v[1] = v[1] / len;

v[2] = v[2] / len;}

//Получить матрицу поворота вокруг вектора v на угол angle

void matrix_rotation(double v[3], double angle, double M[3][3]){

double cosa, sina;

double x, y, z;

cosa = cos(angle*DTR);

sina = sin(angle*DTR);

x = v[0];

y = v[1];

z = v[2];

M[0][0] = cosa + (1 - cosa)*x*x;

M[0][1] = (1 - cosa)*x*y - sina*z;

M[0][2] = (1 - cosa)*x*z + sina*y;

M[1][0] = (1-cosa)*y*x + sina*z;

M[1][1] = cosa + (1 -cosa)*y*y;

M[1][2] = (1 - cosa)*y*z - sina*x;

M[2][0] = (1-cosa)*z*x - sina*y;

M[2][1] = (1 - cosa)*z*y + sina*x;

M[2][2] = cosa+(1 - cosa)*z*z;}

//Расчет угла места

// OP - векртор в J2000 из центра Луны в точку на поверхности Луны;

// OK - вектор в J2000 из центра Луны в координаты КА на гало-орбите.

double elevation_angle(double * OP, double * OK)

{int k;

double lenOP, lenOK, lenPK, betta, gamma, alpha;

double PK[3];// вектор от точки на поверхности луны в положение КА

lenOP = length(OP);//длина вектора OP

lenOK = length(OK);//длина вектора OK

for (k = 0; k < 3; k++) PK[k] = OK[k] - OP[k];

lenPK = length(PK);//длина вектора PK

betta = angle(OP, OK);//угол между векторами OP и OK

gamma = angle(PK, OK);//угол между векторами PK и OK

alpha = 90.e0 - gamma - betta;

return alpha;}

//умножение вектора на матрицу

void vec_on_matrix(double m1[3], double m2[3][3], double m3[3]){

m3[0] = m1[0] * m2[0][0] + m1[1] * m2[0][1] + m1[2] * m2[0][2];

m3[1] = m1[0] * m2[1][0] + m1[1] * m2[1][1] + m1[2] * m2[1][2];

m3[2] = m1[0] * m2[2][0] + m1[1] * m2[2][1] + m1[2] * m2[2][2];

}

//Траспонирование матрицы

void trans_matrix(double M[3][3]){

double M2[3][3];

int i, j;

for (i = 0; i < 3; i++)

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

M2[i][j] = M[j][i];

for (i = 0; i < 3; i++)

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

M[i][j] = M2[i][j];}

void martix_J_to_SG(double Rj[3], double M[3][3], int time[6]){

double dt, T, alph, bet, w, e[13];

int p = 10;

double gaml = 13.17635815;

dt = DateTimeToJulDate(time[0], time[1], time[2], time[3], time[4], time[5]) - 2415020.0;

e[0] = (125.045 - 0.0529921*dt)*DTR;

e[1] = (250.089 - 0.1059842*dt)*DTR;

e[2] = (260.008 + 13.0120009*dt)*DTR;

e[3] = (176.625 + 13.3407154*dt)*DTR;

e[4] = (357.529 + 0.9856003*dt)*DTR;

e[5] = (311.589 + 26.4057084*dt)*DTR;

e[6] = (134.963 + 13.0649930*dt)*DTR;

e[7] = (276.617 + 0.3287146*dt)*DTR;

e[8] = (34.226 + 1.7484877*dt)*DTR;

e[9] = (15.134 - 0.1589763*dt)*DTR;

e[10] = (119.743 + 0.0036096*dt)*DTR;

e[11] = (239.961 + 0.1643573*dt)*DTR;

e[12] = (25.053 + 12.9590088*dt)*DTR;

T = dt / 36525;

alph = (269.9949 + 0.0031*T - 3.8787*sin(e[0]) - 0.1204*sin(e[1]) + 0.07*sin(e[2]) - 0.0172*sin(e[3]) + 0.0072*sin(e[5]) - 0.0052*sin(e[9]) - 0.0043*sin(e[12]))*DTR;

bet = (66.5392 + 0.013*T + 1.5419*cos(e[0]) + 0.0239*cos(e[1]) - 0.0278*cos(e[2]) + 0.0068*cos(e[3]) - 0.0029*cos(e[5]) + 0.0009*cos(e[6]) + 0.0008*cos(e[9]) - 0.0009*cos(e[12]))*DTR;

w = (38.3213 + gaml*dt + gaml + dt - 1.4*(p, (-12))*pow(dt, 2) + 3.561*sin(e[0]) + 0.1208*sin(e[1]) - 0.0642*sin(e[2]) + 0.0158*sin(e[3]) + 0.0252*sin(e[4]) - 0.0066*sin(e[5]) - 0.0047*sin(e[6]) - 0.0046*sin(e[7]) + 0.0028*sin(e[8]) + 0.0052*sin(e[9]) + 0.004*sin(e[10]) + 0.0019*sin(e[11]) - 0.0044*sin(e[12]))*DTR;

M[0][0] = -sin(alph)*cos(w) - cos(alph)*sin(bet)*sin(w);

M[0][1] = cos(alph)*cos(w) - sin(alph)*sin(bet)*sin(w);

M[0][2] = sin(w)*cos(bet);

M[1][0] = sin(alph)*sin(w) - cos(alph)*sin(bet)*cos(w);

M[1][1] = -cos(alph)*sin(w) - sin(alph)*sin(bet)*cos(w);

M[1][2] = cos(w)*cos(bet);

M[2][0] = cos(alph)*cos(bet);

M[2][1] = sin(alph)*cos(bet);

M[2][2] = sin(bet);}

/*Перевод из СК SG(декартовых) в сферические*/

void SG_to_S(double dec[3], double sph[3])

{double x, y, z;

x = dec[0];

y = dec[1];

z = dec[2];

sph[0] = RTD*atan2(x, y); // lat

sph[1] = RTD*atan2(sqrt(x*x + y*y),z); // lon

sph[2] = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); }

void J_to_S(double Rj[3], double Rs[3], int time[6]){

double Rsg[3];

J_to_SG(Rj, Rsg, time);

SG_to_S(Rsg, Rs);}

void S_to_J(double Rs[3], double Rj[3], int time[6]){

double Rsg[3];

S_to_SG(Rs, Rsg);

SG_to_J(Rsg, Rj, time);}

int read_line(FILE * STREAM, char *line){

int i;

if (fgets(line, 256, STREAM) == NULL)

return (-1);

i = strlen(line);

while (i>0 && (line[i - 1] == '\r' || line[i - 1] == '\n')){

i--; line[i] = '\0';}

return i;}

Файл «moon.cpp».

#define _CRT_SECURE_NO_DEPRECATE

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include "moon.h"

int moon_height_read(char*file,

double lat_min, double lat_max,

double lon_min, double lon_max,

int lat_N, int lon_N, moon_height**_Height){

FILE * fd;//o-файл, откуда читаем

int i, j, i2,j2;

double h;

moon_height*Height;

Height = (moon_height*)calloc(1, sizeof(moon_height));

Height->lat_min = lat_min;

Height->lat_max = lat_max;

Height->lon_min = lon_min;

Height->lon_max = lon_max;

Height->lat_N = abs(lat_N);

Height->lon_N = abs(lon_N);

Height->height = (double**)calloc(Height->lat_N , sizeof(double*));

for (i = 0; i < Height->lat_N;i++)

Height->height[i] = (double*)calloc(Height->lon_N, sizeof(double));

if ((fd = fopen(file, "rt")) == 0){

printf("Can't read moon file '%s'\n", file);

return 0;}

for (i = 0; i < Height->lat_N; i++){

for (j = 0; j < Height->lon_N; j++){

i2 = (lat_N < 0) ? (Height->lat_N - 1 - i) : i;

j2 = (lon_N < 0) ? (Height->lon_N - 1 - j) : j;

fscanf(fd, "%lf", &h);

Height->height[i2][j2] = h;}}

*_Height = Height;

printf("Readed file moon. Read %dx%d from lat[%lf, %lf] x lon[%lf, %lf].\n", Height->lat_N,Height->lon_N,

Height->lat_min,Height->lat_max,Height->lon_min,Height->lon_max);

fclose(fd);

return 0;}

double moon_height_get(moon_height*Height, double lat, double lon){

int iLat, iLon;

iLat = (int)(Height->lat_N * (lat - Height->lat_min) / (Height->lat_max - Height->lat_min));

if (iLat < 0) iLat = 0;

if (iLat >= Height->lat_N) iLat = Height->lat_N - 1;

iLon = (int)(Height->lon_N * (lon - Height->lon_min) / (Height->lon_max - Height->lon_min));

if (iLon < 0) iLon = 0;

if (iLon >= Height->lon_N) iLon = Height->lon_N - 1;

return Height->height[iLat][iLon];}

Файл «orbita.cpp».

#define _CRT_SECURE_NO_DEPRECATE

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include "orbita.h"

// Считать орбиту спутника

int orbit_read(char*file, orbita_point**_KA,int*_NKA){

FILE * o;//o-файл, откуда читаем

orbita_point*KA;

int NKA;

int NKA_allocated;

int sch = 0, Nstr = 0;// sch-счетчик для case; Nstr-номер строки

char string[301], *pch;

if ((o = fopen(file, "rt")) == 0){

printf("Can't read orbit file '%s'\n", file);

return 0;}

KA = 0;

NKA = 0;

NKA_allocated = 0;

//считываем координаты KA в СК J2000 и время из "ОРБИТА.txt"

// пропускаем заголовки

fgets(string, 300, o);

while (!feof(o)){

fgets(string, 300, o);

if (string != NULL){

if (NKA + 1 >= NKA_allocated){

if (NKA_allocated == 0){

NKA_allocated = 16;

KA = (orbita_point*)calloc(NKA_allocated, sizeof(orbita_point));}

else{

while (NKA + 1 >= NKA_allocated) NKA_allocated = NKA_allocated * 2;

KA = (orbita_point*)realloc(KA, NKA_allocated * sizeof(orbita_point));}}

pch = strtok(string, " ");

KA[NKA].coords[0] = atof(pch); // X

while (pch != NULL){

if (sch<5) pch = strtok(NULL, " ");

if ((sch>5) && (sch < 9)) pch = strtok(NULL, " .");

if (sch >= 9) pch = strtok(NULL, ":");

switch (sch) {

case 0: KA[NKA].coords[1] = atof(pch); break; // Y

case 1: KA[NKA].coords[2] = atof(pch); break; // Z

case 6: KA[NKA].time[2] = atoi(pch); break; // Day

case 7: KA[NKA].time[1] = atoi(pch); break; // Month

case 8: KA[NKA].time[0] = atoi(pch); break; // Year

case 9: KA[NKA].time[3] = atoi(pch); break; // Hour

case 10: KA[NKA].time[4] = atoi(pch); break; // Min

case 11: KA[NKA].time[5] = atoi(pch); break; // Sec

}

sch++;}

sch = 0;

NKA++;

Nstr++;}}

printf("Readed file orbit. Count of points %d.\n", NKA);

*_KA = KA;

*_NKA = NKA;

return 0;}

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