Skip to content

Instantly share code, notes, and snippets.

@hdf

hdf/RKF7_4D_heterokl_Dll.cpp

Last active July 29, 2018 14:08
Show Gist options
  • Save hdf/2a4f01f24c1f74bd48252ee297ed53b2 to your computer and use it in GitHub Desktop.
Save hdf/2a4f01f24c1f74bd48252ee297ed53b2 to your computer and use it in GitHub Desktop.
Download ZIP
Judit1
Raw
RKF7_4D_heterokl_Dll.cpp
// Ennek a fájlnak a letöltési linkje:
// https://gist.github.com/hdf/2a4f01f24c1f74bd48252ee297ed53b2
// Felfedezendő terület vizuális paraméter generáló: https://gist.github.com/hdf/63f5f088d11039972778
// Linux alatti építés és futtatás (ha nincs fönt g++: "sudo apt-get install g++"):
// g++ -O3 -fopenmp RKF7_4D_heterokl_Dll.cpp -o judit1
// chmod +x ./judit1
// ./judit1 -MMM 10 -t 1300 -o eredmeny.txt
// less ./eredmeny.txt
// (less -ből való kilépéshez nyomd meg a q billentyűt.)
// DLL -é (vagy .so -vá) fordítás g++ -al:
// g++ -O3 -fopenmp -fpic -shared -static RKF7_4D_heterokl_Dll.cpp -o judit1.dll
// Linuxon:
// g++-5 -O3 -fopenmp -fpic -shared RKF7_4D_heterokl_Dll.cpp -o judit1.so
// Teszteléshez: https://gist.github.com/hdf/116176f0de96b87cc557
// Adatbázis rendezése Linuxon:
// sort -t $'\t' -k 17 -g tajkep_.txt -o tajkep_sorted.txt
// Windowson ha van fent mingw64 msys2, .bat fájlba ezt:
// @echo off
// setlocal ENABLEDELAYEDEXPANSION
// set dir=%cd:\=/%
// set dir=%dir::=%
// sh -login -c "cd /%dir% && sort -t $'\t' -k 17 -g tajkep_.txt -o tajkep_sorted.txt"
//A Föld-Hold-műhold probléma "teljesen" valós modellje {változó időléptetésű Runge-Kutta 7(8)-cal, egy kezdőpontból}
//A Hold és a Nap koordínátáit is a mozgásegyenletekből kapjuk
//geocentrikus derékszőgű ekliptikai koordinátarendszerben vagyunk, síkbeli modell, a z koordinátát ehanyagoljuk
//nemcsak a múhold, hanem a nagy Hold és a Nap pályáját is kirajzoljuk, és az L4, L5 Lagrange pontokét is
//a Nap is, a Hold is, a Föld is és a műhold is kiterjedés nélküli tömegpontnak van tekintve
//Poincaré metszet
//síkbeli modell
//kezelve a szoros megközelítés és az ablak bezárása kilépéskor
//i,j-t ír ki futás közben
#ifdef _MSC_VER
#define _CRT_SECURE_NO_WARNINGS
#endif
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
//#include <string>
#include <time.h>
//#undef _OPENMP
#ifdef _OPENMP
#include <omp.h>
#endif
#ifndef _WIN32
#define __stdcall
#endif
const int STR_MAX = 2048;
const double intervall = 100;
//itt következik a kiindulási tartomány, a fixpont koordinátái: végtelen
//main definitions
const int nmax = 12; //A fázistér dimenziója
const double
//ezek a JPL-ből vett adatok
GE = 2.97521695276e15, //köbkilométer/négyzetnap mértékegységben, a Föld gravitációs állandója*tömege
GS = 9.90693056237e20, //a Nap gravitációs állandója*tömege
GM = 0.0366024999763e15, //G a gravitációs állandó, M a Hold tömege, E a Föld tömege, S a Nap tömege
G = 4.98060960352e-10, //köbkm/négyzetnap, 1 kg tömegűnek véve az űrhajót, az űrhajó gravitációs állandója *1 kg
r0 = 4e5, //= FH km,(közepes Föld-Hold távolság, ezzel dimenziótlanítottam az egyenleteket)
//AU=1.4959787066e8 km
kk = 373.99467665, //átszámítási szorzó tényező AU és FH között
//1 AU = 373.99467665 * FH távolság(közepes Föld-Hold távolság
//A)epocha: A következő kezdeti adatok a JDepoch = 2455999.5, 2012. március 13. 00:00 UT - ra vonatkoznak!!!!
//ezek az űrhajó kezdő koordinátái
x0 = -0.9141820107443692, //a valós L4 pont: -0.9141820107443692, L5 kezdőpontja: 0.3975653477547422
yy0 = 0.0687343088982397, //a valós L4 kezdőpontja: 0.0687343088982397, L5 kezdőpontja: -0.8260719994364822
vx0 = -0.0252733218674244, //az L4 pont kezdősebessége: -0.0252733218674244, L5: 0.2106560466347354 //mértékegység: 1/nap
vy0 = -0.2286530912785001, //az L4 pont kezdősebessége: -0.2286530912785001, L5: 0.0924392068640397 1/nap
//a Hold és a Nap kezdőadatai AU-ban vannak megadva, a műholdé és a Lagrange pontoké r0-lal(=FH)dimenziótlanítottak
//Hold: hx0 = -0.516616663 FH, hyy0 = -0.75733769054 FH, hvx0 = 0.185382725, hvy0 = -0.13621388441 FH/nap
hx0 = -1.381347637397253e-03*kk, //a Hold kezdő koordinátái. a mozgásokat a Földéhez viszonyítjuk, AU-ban
hvx0 = 4.956827900007784e-04*kk,
hyy0 = -2.024995909893635e-03*kk,
hvy0 = -3.642134311498105e-04*kk,
nx0 = 9.857531758686030E-01*kk, //a Nap kezdő koordinátái, geocentrikus ekliptikus koord. rendszer
nvx0 = 2.488499979373501E-03*kk,
nyy0 = -1.272767764913456E-01*kk,
nvy0 = 1.712320940981608E-02*kk,
epszilon = 1e-16;
static void derivalt(/*double t,*/ double *x, double *xp) {
double rFH, rFN, rHN, rF4, rH4, rN4, Fx, Fy;
//Az x nmax elemû vektor által kijelölt helyen a t idõben kiszámítja
//az x deriváltját, és xp-ben visszaadja.
rFH = sqrt(x[5] * x[5] + x[7] * x[7]);
rFN = sqrt(x[9] * x[9] + x[11] * x[11]);
rHN = sqrt((x[9] - x[5])*(x[9] - x[5]) + (x[11] - x[7])*(x[11] - x[7]));
rF4 = sqrt(x[1] * x[1] + x[3] * x[3]);
rH4 = sqrt((x[5] - x[1])*(x[5] - x[1]) + (x[7] - x[3])*(x[7] - x[3]));
rN4 = sqrt((x[9] - x[1])*(x[9] - x[1]) + (x[11] - x[3])*(x[11] - x[3]));
Fx = (GM*x[5] / (rFH*rFH*rFH) + GS*x[9] / (rFN*rFN*rFN) + G*x[1] / (rF4*rF4*rF4)) / (r0*r0*r0);
Fy = (GM*x[7] / (rFH*rFH*rFH) + GS*x[11] / (rFN*rFN*rFN) + G*x[3] / (rF4*rF4*rF4)) / (r0*r0*r0);
xp[1] = x[2];
xp[2] = -(GE*x[1] / (rF4*rF4*rF4) + GM *(x[1] - x[5]) / (rH4*rH4*rH4) +
GS*(x[1] - x[9]) / (rN4*rN4*rN4)) / (r0*r0*r0) - Fx;
xp[3] = x[4];
xp[4] = -(GE*x[3] / (rF4*rF4*rF4) + GM *(x[3] - x[7]) / (rH4*rH4*rH4) +
GS*(x[3] - x[11]) / (rN4*rN4*rN4)) / (r0*r0*r0) - Fy;
xp[5] = x[6];
xp[6] = -(GE*x[5] / (rFH*rFH*rFH) + GS *(x[5] - x[9]) / (rHN*rHN*rHN) +
G*(x[5] - x[1]) / (rH4*rH4*rH4)) / (r0*r0*r0) - Fx;
xp[7] = x[8];
xp[8] = -(GE*x[7] / (rFH*rFH*rFH) + GS *(x[7] - x[11]) / (rHN*rHN*rHN) +
G*(x[7] - x[3]) / (rH4*rH4*rH4)) / (r0*r0*r0) - Fy;
xp[9] = x[10];
xp[10] = -(GE*x[9] / (rFN*rFN*rFN) + GM *(x[9] - x[5]) / (rHN*rHN*rHN) +
G*(x[9] - x[1]) / (rN4*rN4*rN4)) / (r0*r0*r0) - Fx;
xp[11] = x[12];
xp[12] = -(GE*x[11] / (rFN*rFN*rFN) + GM *(x[11] - x[7]) / (rHN*rHN*rHN) +
G*(x[11] - x[3]) / (rN4*rN4*rN4)) / (r0*r0*r0) - Fy;
}
//ez a változó lépésközû Runge Kutta Fehler 7(8)-adrendû módszer
static void RungeKuttaFehlberg7(double *x, /*double t,*/ double dt, double *TTmax) {
//Az x nmax elem– x vektort dt idővel továbbfejleszti az RKF7 módszer segítségével
//A derivalt rutin számítja az x deriváltját.
double k1[nmax + 1], k2[nmax + 1], k3[nmax + 1], k4[nmax + 1], k5[nmax + 1], k6[nmax + 1], k7[nmax + 1], k8[nmax + 1], k9[nmax + 1], k10[nmax + 1], k11[nmax + 1], k12[nmax + 1], k13[nmax + 1], y[nmax + 1];
int i;
derivalt(/*t,*/ x, k1);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt * 2 * k1[i] / 27;
derivalt(/*t + dt * 2 / 27,*/ y, k2);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*k1[i] / 36 + dt*k2[i] / 12;
derivalt(/*t + dt / 9,*/ y, k3);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(k1[i] / 24 + k3[i] / 8);
derivalt(/*t + dt / 6,*/ y, k4);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(k1[i] * 5 / 12 - k3[i] * 25 / 16 + k4[i] * 25 / 16);
derivalt(/*t + dt * 5 / 12,*/ y, k5);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(k1[i] / 20 + k4[i] / 4 + k5[i] / 5);
derivalt(/*t + dt / 2,*/ y, k6);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(-k1[i] * 25 / 108 + k4[i] * 125 / 108 - k5[i] * 65 / 27 + k6[i] * 125 / 54);
derivalt(/*t + 5 * dt / 6,*/ y, k7);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(k1[i] * 31 / 300 + k5[i] * 61 / 225 - k6[i] * 2 / 9 + k7[i] * 13 / 900);
derivalt(/*t + dt / 6,*/ y, k8);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(k1[i] * 2 - k4[i] * 53 / 6 + k5[i] * 704 / 45 - k6[i] * 107 / 9 + k7[i] * 67 / 90 + k8[i] * 3);
derivalt(/*t + dt * 2 / 3,*/ y, k9);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(-k1[i] * 91 / 108 + k4[i] * 23 / 108 - k5[i] * 976 / 135 + k6[i] * 311 / 54 - k7[i] * 19 / 60 + k8[i] * 17 / 6 - k9[i] / 12);
derivalt(/*t + dt / 3,*/ y, k10);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(k1[i] * 2383 / 4100 - k4[i] * 341 / 164 + k5[i] * 4496 / 1025 - k6[i] * 301 / 82 + k7[i] * 2133 / 4100 + k8[i] * 45 / 82 + k9[i] * 45 / 164 + k10[i] * 18 / 41);
derivalt(/*t + dt,*/ y, k11);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(k1[i] * 3 / 205 - k6[i] * 6 / 41 - k7[i] * 3 / 205 - k8[i] * 3 / 41 + k9[i] * 3 / 41 + k10[i] * 6 / 41);
derivalt(/*t,*/ y, k12);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(-k1[i] * 1777 / 4100 - k4[i] * 341 / 164 + k5[i] * 4496 / 1025 - k6[i] * 289 / 82 + k7[i] * 2193 / 4100 + k8[i] * 51 / 82 + k9[i] * 33 / 164 + k10[i] * 12 / 41 + k12[i]);
derivalt(/*t + dt,*/ y, k13);
for (i = 1; i <= nmax; i++)
x[i] = x[i] + dt*(k6[i] * 34 / 105 + k7[i] * 9 / 35 + k8[i] * 9 / 35 + k9[i] * 9 / 280 + k10[i] * 9 / 280 + k12[i] * 41 / 840 + k13[i] * 41 / 840);
*TTmax = 0;
double TT;
for (i = 1; i <= nmax; i++) {
TT = fabs(k1[i] + k11[i] - k12[i] - k13[i])*dt * 41 / 840;
if (TT > *TTmax)
*TTmax = TT;
}
}
static double Nyereghalmaz(long bei, long bej, long n, char *sout, double *x, size_t xs, double days, void (__stdcall *mgdFunc)(const char *msg, char *str), bool (__stdcall *mgdFuncAbort)()) {
double t = 0, dt, TTmax, next_t = intervall;
char s[STR_MAX] = "";
sprintf(s, "%18.16f \t %18.16f \t %18.16f \t %18.16f \t", x[1], x[2], x[3], x[4]); //innét indul a trajektória t=0-ban
do {
//innét van a nyereghalmazos beszúrás, az days napig bent maradó trajektóriák kezdőpontját, egy közbenső időpontbeli pontját és végpontját írja ki
if (mgdFunc != NULL) {
char s2[STR_MAX] = "";
#ifdef _OPENMP
sprintf(s2, "%s %8.6f \t %18.16f \t %18.16f \t %18.16f \t %18.16f \t %ld \t %ld \t %ld \t %d", s, t, x[1], x[2], x[3], x[4], bei, bej, n, omp_get_thread_num());
#pragma omp critical
mgdFunc("partial", s2);
#else
sprintf(s2, "%s %8.6f \t %18.16f \t %18.16f \t %18.16f \t %18.16f \t %ld \t %ld \t %ld \t %d", s, t, x[1], x[2], x[3], x[4], bei, bej, n, 0);
mgdFunc("partial", s2);
#endif
}
if (t >= days) {
#ifdef _OPENMP
sprintf(s + strlen(s), " %8.6f \t %18.16f \t %18.16f \t %18.16f \t %18.16f \t %ld \t %ld \t %ld \t %d", t, x[1], x[2], x[3], x[4], bei, bej, n, omp_get_thread_num());
#else
sprintf(s + strlen(s), " %8.6f \t %18.16f \t %18.16f \t %18.16f \t %18.16f \t %ld \t %ld \t %ld \t %d", t, x[1], x[2], x[3], x[4], bei, bej, n, 0);
#endif
sprintf(sout, "%s\n", s);
break;
} else if (t >= next_t) { // megadott időközönként hol tart a trajektória (nyereghalmaz)
sprintf(s + strlen(s), " %8.6f \t %18.16f \t %18.16f \t %18.16f \t %18.16f \t", t, x[1], x[2], x[3], x[4]);
next_t = t + intervall;
}
//idáig van a nyereghalmazos beszúrás
dt = 0.01;
do {
double x1[nmax + 1];
memcpy(x1, x, xs);
RungeKuttaFehlberg7(x1, /*t,*/ dt, &TTmax);
//ha belép az alábbi tartományba, akkor nem folytatjuk, azaz, ha 1.5 sugárnál távolabb kerüla Földtől az űrhajó
if (sqrt(x[1] * x[1] + x[3] * x[3]) > 1.5)
return t;
//ha túl közel kerül a holdhoz az űrhajó, álljon le a program
if (sqrt((x[1] - x[5])*(x[1] - x[5]) + (x[3] - x[7])*(x[3] - x[7])) < 0.000025) {
sprintf(sout, "%s\n", s);
break;
}
if (TTmax > epszilon)
dt = 0.9*dt*pow(epszilon / TTmax, 1 / 8);
else
memcpy(x, x1, xs);
} while (TTmax > epszilon);
t += dt;
} while (mgdFuncAbort == NULL || !mgdFuncAbort());
return -1; //ahol -1-gyel lép ki a program, nem lép ki a trajektória a megadott idõn belül
}
extern "C" { // Ez azért kell, hogy ezek a függvények elérhetőek legyenek JS -ben
void calculate(long i, long j, long n, long MMM, long NNN, double days, char *out, double xstart, double ystart, double xres, double yres, void (__stdcall *mgdFunc)(const char *msg, char *str), bool (__stdcall *mgdFuncAbort)()) {
double x[nmax + 1]; //rFH a Nap-Jupiter távolsága, stb., a 4-es index a tömeg nélküli negyedik testet jelenti
x[1] = xstart + xres * i / MMM;
x[3] = ystart + yres * j / NNN;
x[2] = vx0;
x[4] = vy0;
x[5] = hx0;
x[6] = hvx0;
x[7] = hyy0;
x[8] = hvy0;
x[9] = nx0;
x[10] = nvx0;
x[11] = nyy0;
x[12] = nvy0;
Nyereghalmaz(i, j, n, out, x, sizeof(x), days, mgdFunc, mgdFuncAbort);
}
#if defined(_WIN32) && defined(_DLL)
__declspec(dllexport)
#endif
void asyncCalculator(void (__stdcall *mgdFunc)(const char *msg, char *str), bool (__stdcall *mgdFuncAbort)(), bool partial, long begin, long end, long MMM, long NNN, double days, double xstart, double ystart, double xres, double yres) {
if (end < 1 || end > ((MMM + 1)*(NNN + 1)))
end = ((MMM + 1)*(NNN + 1));
clock_t start = clock();
#ifdef _OPENMP
#pragma omp parallel default(none) firstprivate(MMM, NNN, days, xstart, ystart, xres, yres, begin, end, partial) shared(mgdFunc, mgdFuncAbort)
{
#pragma omp for schedule(dynamic)
#endif
//itt következik a kiindulási tartomány,
//for(long i = 0; i <= MMM; i++) {
//for(long j = 0; j <= NNN; j++) {
for (long n = begin; n < end; n++) { // 1 dimenzióssá lapított 2 dimenziós tömb, mivel OpenMP 2.0 nem tud több dimenziós tömböt párhuzamosítani
long i = n / (NNN + 1);
long j = n % (NNN + 1);
char s[STR_MAX] = "";
calculate(i, j, n, MMM, NNN, days, s, xstart, ystart, xres, yres, ((partial)?mgdFunc:NULL), mgdFuncAbort);
#ifdef _OPENMP
#pragma omp critical
if (s[0] != '\0') mgdFunc("stable", s);
#endif
//printf("Nyereghalmaz meghívva: %ld\t%ld\tSzál: %d\n", i, j, omp_get_thread_num());
//sprintf(out + strlen(out), "%s", s);
//}
}
#ifdef _OPENMP
}
#endif
float elapsed = ((float)(clock() - start)) / CLOCKS_PER_SEC;
#if defined(_OPENMP) && defined(_WIN32)
//elapsed = elapsed / (float)omp_get_num_procs(); // Not needed when called from C#???
#endif
char s[16] = "";
sprintf(s, "%.3f", elapsed);
mgdFunc("finished", s);
}
#if defined(_WIN32) && defined(_DLL)
__declspec(dllexport)
#endif
float calculator(char *out, bool dll, long begin, long end, long MMM, long NNN, double days, double xstart, double ystart, double xres, double yres) {
if (end < 1 || end > ((MMM + 1)*(NNN + 1)))
end = ((MMM + 1)*(NNN + 1));
FILE *stream = NULL;
if (!dll) stream = fopen(out, ((begin > 0)?"a":"w"));
clock_t start = clock();
#ifdef _OPENMP
//std::string results = "";
#pragma omp parallel default(none) firstprivate(MMM, NNN, days, xstart, ystart, xres, yres, begin, end, dll) shared(stream, out)
{
#pragma omp for schedule(dynamic)
#endif
//itt következik a kiindulási tartomány,
//for(long i = 0; i <= MMM; i++) {
//for(long j = 0; j <= NNN; j++) {
for (long n = begin; n < end; n++) { // 1 dimenzióssá lapított 2 dimenziós tömb, mivel OpenMP 2.0 nem tud több dimenziós tömböt párhuzamosítani
long i = n / (NNN + 1);
long j = n % (NNN + 1);
char s[STR_MAX] = "";
calculate(i, j, n, MMM, NNN, days, s, xstart, ystart, xres, yres, NULL, NULL);
if (!dll) {
#ifdef _OPENMP
#pragma omp critical
//results += s; // Ez egyszerre csak 1 szálon történhet
printf("Nyereghalmaz meghívva: %ld\t%ld\tSzál: %d\n", i, j, omp_get_thread_num());
#else
printf("Nyereghalmaz meghívva: %ld\t%ld\n", i, j);
#endif
fprintf(stream, "%s", s);
} else sprintf(out + strlen(out), "%s", s);
//}
}
#ifdef _OPENMP
}
#endif
float elapsed = ((float)(clock() - start)) / CLOCKS_PER_SEC;
#ifdef _OPENMP
#ifndef _WIN32
elapsed = elapsed / (float)omp_get_num_procs();
#endif
//fprintf(stream, "%s", results.c_str());
#endif
if (!dll) fclose(stream);
return elapsed;
}
int main(int argc, char *argv[]) {
#ifdef _WIN32
system("chcp 65001"); // Ez azért kell, hogy windowson is működjön az ékezetes betű kiírás
#endif
char fileName[0x100] = "tajkep_.txt";
long MMM = 10;
long NNN = 10;
double days = 200;
double xstart = -1.15;
double ystart = -0.4;
double xres = 0.45;
double yres = 0.65;
long begin = 0;
long end = 0;
// Kapcsolók feldolgozása
for (int i = 1; i < argc; i += 2) {
if (strcmp(argv[i], "-o") == 0)
sprintf(fileName, "%s", argv[i + 1]);
else if (strcmp(argv[i], "-MMM") == 0)
MMM = strtol(argv[i + 1], NULL, 10);
else if (strcmp(argv[i], "-NNN") == 0)
NNN = strtol(argv[i + 1], NULL, 10);
else if (strcmp(argv[i], "-t") == 0)
days = strtod(argv[i + 1], NULL); //ennyi ideig tart a futás, a mértékegység: [nap]
else if (strcmp(argv[i], "-xs") == 0)
xstart = strtod(argv[i + 1], NULL);
else if (strcmp(argv[i], "-ys") == 0)
ystart = strtod(argv[i + 1], NULL);
else if (strcmp(argv[i], "-aw") == 0)
xres = strtod(argv[i + 1], NULL);
else if (strcmp(argv[i], "-ah") == 0)
yres = strtod(argv[i + 1], NULL);
else if (strcmp(argv[i], "-s") == 0)
begin = strtol(argv[i + 1], NULL, 10);
else if (strcmp(argv[i], "-e") == 0)
end = strtol(argv[i + 1], NULL, 10);
else {
printf("Usage: %s [options]\nOptions:\n"
" -MMM matrix_y_size\t\t(default: %ld)\n"
" -NNN matrix_x_size\t\t(default: %ld)\n"
" -o output_file\t\t(default: \"%s\")\n"
" -t days_to_run\t\t(default: %.3f)\n"
" -xs area_x_offset\t\t(default: %.3f)\n"
" -ys area_y_offset\t\t(default: %.3f)\n"
" -aw area_width\t\t(default: %.3f)\n"
" -ah area_height\t\t(default: %.3f)\n"
" -s begin_at_particle\t\t(default: %ld)\n"
" -e end_after_particle\t\t(default: %ld)\n",
argv[0], MMM, NNN, fileName, days, xstart, ystart, xres, yres, begin, end);
return 0;
}
}
printf("Futott: %.3f másodpercig\n",
calculator(fileName, false, begin, end, MMM, NNN, days, xstart, ystart, xres, yres));
//getchar();
return 0;
}
}
Raw
RKF7_terbeli_bari_Fold_Hold_heterokl_6D.cpp
//forrásfájl: RKF7_4D_bari_Naprendszer_kontroll_param_T.cpp
//A epocha
//térbeli modell, a sebességeket is léptetjük, azért 6D!!!
//az új, a Föld-Hold-űrhajó tkp-jával is jó modell (azaz a Nap nélkül)
//űrhajó, Nap, Föld, Hold modell, ezek közös tkp-hoz képest vannak a koordináták
//a kezdeti feltételekből ki kell vonni a tkpx0-akat, számolás így folyik. majd a végén ez eredményhez hozzá kell adni a tkpx0-akat. Így lesz konzisztens a kimenő adat a bemenőekkel.
//az aktualis tömegek baricentruma a középpont, ekliptikai koordinátarendszerben vagyunk,
//a kezdokoordinatak a JPL-bol vannak veve, a Naprendszer baricentrumahoz kepest. Ezek a koordinatak at vannak szamolva az aktualis tomegek barucentrumahoz kepesti ekliptikai koordibatarendszerre
//az egitestek szam tetszolegesen valtoztathato
//valamennyi test kiterjedes nelkuli tomegpontnak van tekintve
//külön fájlokba printeli a különböző időpontokhoz tartozó koordinátákat
// https://gist.github.com/hdf/2a4f01f24c1f74bd48252ee297ed53b2#file-rkf7_terbeli_bari_fold_hold_heterokl_6d-cpp
// g++ -O3 -fopenmp RKF7_terbeli_bari_Fold_Hold_heterokl_6D.cpp -o judit1
#ifdef _MSC_VER
#define _CRT_SECURE_NO_WARNINGS
#endif
#ifdef __unix
#define fopen_s(pFile,filename,mode) ((*(pFile))=fopen((filename),(mode)))==NULL
#endif
// C RunTime Header Files
#include <stdlib.h>
#include <stdio.h>
#include <memory.h>
#include <math.h>
#include <time.h>
#include <omp.h>
#define PI 3.1415926535897932384626433832795
// main definitions
const int nmax = 24; // A fazister dimenzioja}
const int N = 4; // az egitestek szama
//double xMin = -3.5, xMax =3.5 ;// ezek a hatarok xy-ban, átírhatók, ha nem jók
//double yMin = -2, yMax = 2;
//double xMin = -0.00875, xMax = 0.00875; // ezek a hatarok xy-ban, átírhatók, ha nem jók
//double yMin = -0.005, yMax = 0.005;
const double
//ezek a JPL-bol vett adatok
G = 4.98060960352e-10, //kobkm/kg/negyzetnap
GS = 9.90693056237e20, //kobkm/kg/negyzetnap
GJ = 9.45384450672e17, //kobkm/kg/negyzetnap
GSzat = 2.83154948987e17,
GMars = 3.19711546368e14,
GFold = 2.97553634058e15,
mnap = 1.988544e30, //kg
mhold = 0.007349e25,
mfold = 0.597219e25,
mmerkur = 3.302e23,
mvenusz = 48.685e23,
mmars = 6.4185e23,
mjupiter = 1.89813e27,
mszaturnusz = 5.68319e26,
muranusz = 86.8103e24,
mneptunusz = 102.41e24,
ajup = 5.20336301, //AU, a Jupiter nagytengelye
m[N + 1] = { 0, 1, 1.988544e30, 0.597219e25, 0.007349e25 }, // 1.89813e27}, //
//m[N+1] = {0, 1, 1.988544e30, 1898.13e24, 5.68319e26, 0.597219e25, 0.007349e25, 3.302e23, 48.685e23, 6.4185e23, 86.8103e24, 102.41e24},
// a tomegek, kg: urhajo, Nap, Jupiter, Szaturnusz,Fold, Hold, Merkur, Venusz, Mars, Uranusz, Neptunusz sorrendben
//r0 = 4e5, //= FH (kozepes Fold-Hold tavolsag, ezzel dimenziotlanitottam az egyenleteket)
r0 = 1.4959787066e8, //1 AU, ezzel dimenziotlanitottam az egyenleteket
//AU=1.4959787066e8 km
//kk = 373.99467665, //atszamitasi szorzo tenyezo AU es FH kozott
// 1 AU = 373.99467665 * FH tavolsag (kozepes Fold-Hold tavolsag
// A-7) epocha: A kovetkezo kezdeti adatok a 2457976.458333333 = A.D. 2017-Aug-10 23:00:00.0000 TDB -ra vonatkoznak AU-ban, illetve AU/nap-ban
fx0 = 7.586742351597854E-01, // t=0-ban a Föld kezdő koordinátái
fvx0 = 1.117119817220456E-02,
fyy0 = -6.695310414037162E-01,
fvy0 = 1.277008443245662E-02,
fz0 = -1.116185678411723E-04,
fvz0 = 7.329437509563111E-08,
/*
x0 = 0.0014154757948744 + fx0, // az űrhajó valódi kezdő koordinátái a Föld-Hold L4 pontjából,
vx0= -0.0005050483356001 + fvx0,
yy0=0.0021277694610626 + fyy0,
vy0=0.0003094527247440 + fvy0,
z0 =-0.0002362880499492 + fz0,
vz0 =0.0000043113143139 + fvz0,
*/
x0 = 0.0011423780154981 + fx0, //t0=0-ben az űrhajó valódi kezdő koordinátái a Föld-Hold L5 pontjából ,
vx0 = 0.0005190996126544 + fvx0,
yy0 = -0.0022955927738733 + fyy0,
vy0 = 0.0002809997586600 + fvy0,
z0 = 0.0001096947314251 + fz0,
vz0 = -0.0000492760119127 + fvz0,
hx0 = 7.612320690990555E-01, // t=0-ban a Hold kezdő koordinátái
hvx0 = 1.118525059202575E-02,
hyy0 = -6.696991341636114E-01,
hvy0 = 1.336053848828723E-02,
hz0 = -2.382559188457015E-04,
hvz0 = -4.487039784450909E-05,
nx0 = 2.555639022341175E-03, // t=0-ban a Nap kezdő koordinátái
nvx0 = -4.709048117389841E-06,
nyy0 = 5.796156280250261E-06,
nvy0 = 5.791155635556258E-06,
nz0 = -1.378381642936898E-04,
nvz0 = 1.054909391718740E-07,
M = m[1] + m[2] + m[3] + m[4],
tkpx0 = m[1] * x0 / M + m[2] * nx0 / M + m[3] * fx0 / M + m[4] * hx0 / M, // ezek a tomegkozeppont hely- es sebessegkoordinatai
tkpy0 = m[1] * yy0 / M + m[2] * nyy0 / M + m[3] * fyy0 / M + m[4] * hyy0 / M,
tkpz0 = m[1] * z0 / M + m[2] * nz0 / M + m[3] * fz0 / M + m[4] * hz0 / M,
tkpvx0 = m[1] * vx0 / M + m[2] * nvx0 / M + m[3] * fvx0 / M + m[4] * hvx0 / M,
tkpvy0 = m[1] * vy0 / M + m[2] * nvy0 / M + m[3] * fvy0 / M + m[4] * hvy0 / M,
tkpvz0 = m[1] * vz0 / M + m[2] * nvz0 / M + m[3] * fvz0 / M + m[4] * hvz0 / M,
epszilon = 1e-16,
delta = 0.01,
vxmetsz = 0, //0.0005095389099076, //ez a Poincaré metszet magassága, vx-szel metszem az x-y síkot, (vx0-fvx0)
tMax = 8; //3650; //ennyi ideig tart a futás, a mertekegyseg: [nap] 164250 = 450 keringés
///double x[nmax + 1], r[N + 1][N + 1], rN1, rJ1, rNJ, vkep, vxkep, vykep, tkpx, tkpy, tkpvx, tkpvy, fell, vxkepell, vykepell, vkepell, Fx, Fy, vsz, lax, lay;
///double xregi[nmax + 1], xmetsz[nmax + 1], FL5x, FL5vx, FL5y, FL5vy, ML5x, ML5vx, ML5y, ML5vy, FL4x, FL4vx, FL4y, FL4vy, ML4x, ML4vx, ML4y, ML4vy;
//double TTmax, TT, dt;
//------------- functions
void derivalt(/*double t,*/ double* x, double* xp) {
// Az x nmax elemu vektor  altal kijelolt helyen a t idoben kiszamitja}
// az x derivaltjat, es xp-ben visszaadja. }
int n, h;
double r[N + 1][N + 1];
for (n = 1; n <= N; n++) {
xp[6 * n - 4] = 0;
xp[6 * n - 2] = 0;
xp[6 * n] = 0;
for (h = 1; h <= N; h++) {
if (n != h){
r[n][h] = sqrt((x[6 * h - 5] - x[6 * n - 5])*(x[6 * h - 5] - x[6 * n - 5]) + (x[6 * h - 3] - x[6 * n - 3])*(x[6 * h - 3] - x[6 * n - 3]) + (x[6 * h - 1] - x[6 * n - 1])*(x[6 * h - 1] - x[6 * n - 1]));
xp[6 * n - 5] = x[6 * n - 4];
xp[6 * n - 4] = xp[6 * n - 4] - (G*m[h] * (x[6 * n - 5] - x[6 * h - 5]) / (r[n][h] * r[n][h] * r[n][h]) / (r0*r0*r0));
xp[6 * n - 3] = x[6 * n - 2];
xp[6 * n - 2] = xp[6 * n - 2] - (G*m[h] * (x[6 * n - 3] - x[6 * h - 3]) / (r[n][h] * r[n][h] * r[n][h]) / (r0*r0*r0));
xp[6 * n - 1] = x[6 * n];
xp[6 * n] = xp[6 * n] - (G*m[h] * (x[6 * n - 1] - x[6 * h - 1]) / (r[n][h] * r[n][h] * r[n][h]) / (r0*r0*r0));
}
}
}
}
//ez a változó lépésközu Runge Kutta Fehler 7(8)-adrendu módszer
void RungeKuttaFehlberg7(double* x, /*double t,*/ double dt, double* TTmax) {
// Az x nmax elemű x vektort dt idővel tov bbfejleszti az RKF7 módszer segítségével}
// A derivalt rutin sz m!tja az x }
// deriv ltj t. }
double k1[nmax + 1], k2[nmax + 1], k3[nmax + 1], k4[nmax + 1], k5[nmax + 1], k6[nmax + 1], k7[nmax + 1], k8[nmax + 1], k9[nmax + 1], k10[nmax + 1], k11[nmax + 1], k12[nmax + 1], k13[nmax + 1], y[nmax + 1];
int i;
derivalt(/*t,*/ x, k1);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt * 2 * k1[i] / 27;
derivalt(/*t + dt * 2 / 27,*/ y, k2);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*k1[i] / 36 + dt*k2[i] / 12;
derivalt(/*t + dt / 9,*/ y, k3);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(k1[i] / 24 + k3[i] / 8);
derivalt(/*t + dt / 6,*/ y, k4);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(k1[i] * 5 / 12 - k3[i] * 25 / 16 + k4[i] * 25 / 16);
derivalt(/*t + dt * 5 / 12,*/ y, k5);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(k1[i] / 20 + k4[i] / 4 + k5[i] / 5);
derivalt(/*t + dt / 2,*/ y, k6);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(-k1[i] * 25 / 108 + k4[i] * 125 / 108 - k5[i] * 65 / 27 + k6[i] * 125 / 54);
derivalt(/*t + 5 * dt / 6,*/ y, k7);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(k1[i] * 31 / 300 + k5[i] * 61 / 225 - k6[i] * 2 / 9 + k7[i] * 13 / 900);
derivalt(/*t + dt / 6,*/ y, k8);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(k1[i] * 2 - k4[i] * 53 / 6 + k5[i] * 704 / 45 - k6[i] * 107 / 9 + k7[i] * 67 / 90 + k8[i] * 3);
derivalt(/*t + dt * 2 / 3,*/ y, k9);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(-k1[i] * 91 / 108 + k4[i] * 23 / 108 - k5[i] * 976 / 135 + k6[i] * 311 / 54 - k7[i] * 19 / 60 + k8[i] * 17 / 6 - k9[i] / 12);
derivalt(/*t + dt / 3,*/ y, k10);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(k1[i] * 2383 / 4100 - k4[i] * 341 / 164 + k5[i] * 4496 / 1025 - k6[i] * 301 / 82 + k7[i] * 2133 / 4100 + k8[i] * 45 / 82 + k9[i] * 45 / 164 + k10[i] * 18 / 41);
derivalt(/*t + dt,*/ y, k11);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(k1[i] * 3 / 205 - k6[i] * 6 / 41 - k7[i] * 3 / 205 - k8[i] * 3 / 41 + k9[i] * 3 / 41 + k10[i] * 6 / 41);
derivalt(/*t,*/ y, k12);
for (i = 1; i <= nmax; i++)
y[i] = x[i] + dt*(-k1[i] * 1777 / 4100 - k4[i] * 341 / 164 + k5[i] * 4496 / 1025 - k6[i] * 289 / 82 + k7[i] * 2193 / 4100 + k8[i] * 51 / 82 + k9[i] * 33 / 164 + k10[i] * 12 / 41 + k12[i]);
derivalt(/*t + dt,*/ y, k13);
for (i = 1; i <= nmax; i++)
x[i] = x[i] + dt*(k6[i] * 34 / 105 + k7[i] * 9 / 35 + k8[i] * 9 / 35 + k9[i] * 9 / 280 + k10[i] * 9 / 280 + k12[i] * 41 / 840 + k13[i] * 41 / 840);
*TTmax = 0;
double TT;
for (i = 1; i <= nmax; i++) {
TT = fabs(k1[i] + k11[i] - k12[i] - k13[i])*dt * 41 / 840;
if (TT > *TTmax)
*TTmax = TT;
}
}
FILE *stream0, *stream1, *stream2, *stream3, *stream4, *stream5, *stream6; // deklaracio, stream helyett mas nevet is hasznalhatsz persze, ami Neked tetszik
double Nyereghalmaz(int bei, int bej, int bek, int beii, int bejj, int bekk, double* x, size_t xs) {
long szaml;
double t = 0, dt, TTmax;
//fprintf (stream, " %8.6f %18.16f %18.16f \n",t, xmetsz[1],xmetsz[2]);// (-sqrt(3.0)/2)*x[7] + x[5]/2, (sqrt(3.0)/2)*x[5] + x[7]/2);
char s0[1000] = "";
char s1[1000] = "";
char s2[1000] = "";
char s3[1000] = "";
char s4[1000] = "";
char s5[1000] = "";
char s6[1000] = "";
bool megvolt1 = false;
bool megvolt2 = false;
bool megvolt3 = false;
bool megvolt4 = false;
bool megvolt5 = false;
bool megvolt6 = false;
//xmetsz[1] = 0;
//xmetsz[3] = 0;
//a következő 1 sor a pont kezdeti távolságát adja meg a Földtől, az "A" feltételhez kell
double rr0 = sqrt((x[1] - x[13])*(x[1] - x[13]) + (x[3] - x[15])*(x[3] - x[15]));
for (szaml = 0; t <= tMax; szaml++) {
//innét van a nyereghalmazos beszúrás
if (szaml == 0) //innét indul a trajektória t=0-ban, ez nincs a metszeten
sprintf(s0, "%8.6f \t %8.16f \t %8.16f \t %8.16f \t ", t, x[1] - x[13], x[3] - x[15], x[5] - x[17]); //, bei, bej, bek, beii,bejj, bekk);
//fprintf (stream0, " %8.6f \t %8.16f \t %8.16f \t %8.16f \t\n ",t, x[1]-x[13], x[3]-x[15], x[5] - x[17]);
//sprintf (s, "%8.16f \t %8.16f \t", vx0-fvx0, vy0-fvy0);
double zzz = 1;
if ((t > zzz) && !megvolt1) { // félidokor itt van a trajektória
sprintf(s1, "%8.6f \t %8.16f \t %8.16f \t %8.16f \t ", t, x[1] - x[13], x[3] - x[15], x[5] - x[17]);
//fprintf (stream1, " %8.6f \t %8.16f \t %8.16f \t %8.16f \t\n ",t, x[1]-x[13], x[3]-x[15], x[5] - x[17]);
//sprintf (s+ strlen (s), "%8.6f \t %8.16f \t %8.16f \t %8.16f \t %8.16f \t ", t, x[1], x[3], x[1] - x[9],x[3] - x[11]);
megvolt1 = true;
}
double bbb = 2;
if ((t > bbb) && !megvolt2) { // félidokor itt van a trajektória
sprintf(s2, "%8.6f \t %8.16f \t %8.16f \t %8.16f \t ", t, x[1] - x[13], x[3] - x[15], x[5] - x[17]);
//fprintf (stream2, " %8.6f \t %8.16f \t %8.16f \t %8.16f \t\n ",t, x[1]-x[13], x[3]-x[15], x[5] - x[17]);
//sprintf (s+ strlen (s), "%8.6f \t %8.16f \t %8.16f \t %8.16f \t %8.16f \t ", t, x[1], x[3], x[1] - x[9],x[3] - x[11]);
megvolt2 = true;
}
double rrr = 3;
if ((t > rrr) && !megvolt3) { // t=rrr-kor itt van a trajektória
sprintf(s3, "%8.6f \t %8.16f \t %8.16f \t %8.16f \t ", t, x[1] - x[13], x[3] - x[15], x[5] - x[17]);
//fprintf (stream3, " %8.6f \t %8.16f \t %8.16f \t %8.16f \t\n ",t, x[1]-x[13], x[3]-x[15], x[5] - x[17]);
//sprintf (s+ strlen (s), "%8.6f \t %8.16f \t %8.16f \t %8.16f \t %8.16f \t", t, xmetsz[1], xmetsz[3], xmetsz[2], xmetsz[4]);
//sprintf (s+ strlen (s), "%8.6f \t %8.16f \t %8.16f \t %8.16f \t %8.16f \t ", t, x[1], x[3], x[1] - x[9],x[3] - x[11]);
megvolt3 = true;
}
double yyy = 4;
if ((t > yyy) && !megvolt4){ // a végén itt van a trajektória
sprintf(s4, "%8.6f \t %8.16f \t %8.16f \t %8.16f \t ", t, x[1] - x[13], x[3] - x[15], x[5] - x[17]);
//fprintf (stream4, " %8.6f \t %8.16f \t %8.16f \t %8.16f \t\n ",t, x[1]-x[13], x[3]-x[15], x[5] - x[17]);
//sprintf (s+ strlen (s), "%8.6f \t %8.16f \t %8.16f \t %8.16f \t %8.16f \t %d \t %d \t ", t, xmetsz[1], xmetsz[3],xmetsz[2], xmetsz[4], bei, bej);
//sprintf (s+ strlen (s), "%8.6f \t %8.16f \t %8.16f \t %8.16f \t %8.16f \t ", t, x[1], x[3], x[1] - x[9],x[3] - x[11]);
megvolt4 = true;
}
double xxx = 5;
if ((t > xxx) && !megvolt5){ // a végén itt van a trajektória
sprintf(s5, "%8.6f \t %8.16f \t %8.16f \t %8.16f \t ", t, x[1] - x[13], x[3] - x[15], x[5] - x[17]);
//fprintf (stream5, " %8.6f \t %8.16f \t %8.16f \t %8.16f \t\n ", t,x[1]-x[13], x[3]-x[15], x[5] - x[17]);
//sprintf (s+ strlen (s), "%8.6f \t %8.16f \t %8.16f \t %8.16f \t %8.16f \t %d \t %d \t ", t, xmetsz[1], xmetsz[3],xmetsz[2], xmetsz[4], bei, bej);
//sprintf (s+ strlen (s), "%8.6f \t %8.16f \t %8.16f \t %8.16f \t %8.16f \t ", t, x[1], x[3], x[1] - x[9],x[3] - x[11]);
megvolt5 = true;
}
double sss = 6;
if ((t > sss) && !megvolt6){ // a végén itt van a trajektória
sprintf(s6, "%8.6f \t %8.16f \t %8.16f \t %8.16f \t", t, x[1] - x[13], x[3] - x[15], x[5] - x[17]);
//fprintf (stream6, " %8.6f \t %8.16f \t %8.16f \t %8.16f \t\n ",t, x[1]-x[13], x[3]-x[15], x[5] - x[17]);
//sprintf (s+ strlen (s), "%8.6f \t %8.16f \t %8.16f \t %8.16f \t %8.16f \t %d \t %d \t ", t, xmetsz[1], xmetsz[3],xmetsz[2], xmetsz[4], bei, bej);
//sprintf (s+ strlen (s), "%8.6f \t %8.16f \t %8.16f \t %8.16f \t %8.16f \t ", t, x[1], x[3], x[1] - x[9],x[3] - x[11]);
megvolt6 = true;
}
dt = 0.1;
do {
double x1[nmax + 1];
memcpy(x1, x, xs);
RungeKuttaFehlberg7(x1, /*t,*/ dt, &TTmax);
/*
//ez az elvileg helyes, a Nap-Föld L4, L5 pontja a nagyobbik primary-hez képest
FL5vx = -sqrt(x[10] * x[10] + x[12] * x[12]) * sin(PI / 3 - atan(-x[10] / x[12])); // ez az L5 pont kezdosebessége x irányban
FL5vy = -sqrt(x[10] * x[10] + x[12] * x[12]) * cos(PI / 3 - atan(-x[10] / x[12])); // ez az L5 pont kezdosebessége y irányban
FL4vx = sqrt(x[10] * x[10] + x[12] * x[12]) * sin(PI / 3 - atan(x[10] / x[12])); // ez az L4 pont kezdosebessége x irányban
FL4vy = -sqrt(x[10] * x[10] + x[12] * x[12]) * cos(PI / 3 - atan(x[10] / x[12])); // ez az L4 pont kezdosebessége y irányban
FL5x = (sqrt(3.0) / 2)*x[11] + x[9] / 2; // ez az L5 pont kezdo x koordinátája
FL5y = (-sqrt(3.0) / 2)*x[9] + x[11] / 2; // ez az L5 pont kezdo y koordinátája
FL4x = -(sqrt(3.0) / 2)*x[11] + x[9] / 2; // ez az L4 pont kezdo x koordinátája
FL4y = (sqrt(3.0) / 2)*x[9] + x[11] / 2; // ez az L4 pont kezdo y koordinátája
// ez az elvileg helyes, a Nap-Mars L4, L5 pontja a nagyobbik primary-hez képest
ML5vx = -sqrt(x[18] * x[18] + x[20] * x[20]) * sin(PI / 3 - atan(-x[18] / x[20])); // ez az L5 pont kezdosebessége x irányban
ML5vy = -sqrt(x[18] * x[18] + x[20] * x[20]) * cos(PI / 3 - atan(-x[18] / x[20])); // ez az L5 pont kezdosebessége y irányban
ML4vx = sqrt(x[18] * x[18] + x[20] * x[20]) * sin(PI / 3 - atan(x[18] / x[20])); // ez az L4 pont kezdosebessége x irányban
ML4vy = -sqrt(x[18] * x[18] + x[20] * x[20]) * cos(PI / 3 - atan(x[18] / x[20])); // ez az L4 pont kezdosebessége y irányban
ML5x = (sqrt(3.0) / 2)*x[19] + x[17] / 2; // ez az L5 pont kezdo x koordinátája
ML5y = (-sqrt(3.0) / 2)*x[17] + x[19] / 2; // ez az L5 pont kezdo y koordinátája
ML4x = -(sqrt(3.0) / 2)*x[19] + x[17] / 2; // ez az L4 pont kezdo x koordinátája
ML4y = (sqrt(3.0) / 2)*x[17] + x[19] / 2; // ez az L4 pont kezdo y koordinátája
*/
// ha belép az alábbi tartományba, akkor nem folytatjuk
if (sqrt((x[1] - x[19])*(x[1] - x[19]) + (x[3] - x[21])*(x[3] - x[21]) + (x[5] - x[23])*(x[5] - x[23])) < 0.0004 || rr0*1.5 < sqrt((x[1] - x[13])*(x[1] - x[13]) + (x[3] - x[15])*(x[3] - x[15]) + (x[5] - x[17])*(x[5] - x[17])) || sqrt((x[1] - x[13])*(x[1] - x[13]) + (x[3] - x[15])*(x[3] - x[15] + (x[5] - x[17])*(x[5] - x[17]))) < rr0*0.5)
return t;
if (TTmax > epszilon)
dt = 0.9*dt*pow(epszilon / TTmax, 1 / 8);
else
memcpy(x, x1, xs);
} while (TTmax > epszilon);
t += dt;
}
if (megvolt6 == true) {
#pragma omp critical
{
fprintf(stream0, "%s\n", s0);
fprintf(stream1, "%s\n", s1);
fprintf(stream2, "%s\n", s2);
fprintf(stream3, "%s\n", s3);
fprintf(stream4, "%s\n", s4);
fprintf(stream5, "%s\n", s5);
fprintf(stream6, "%s\n", s6);
}
}
return -1; //ahol -1-gyel lép ki a program, nem lép ki a trajektória a megadott idon belül
}
void DrawContent() {
fopen_s(&stream0, "A_proba_0.txt", "w"); // megnyitas irasra, cdata.txt a file neve, ha nem a program mellett van (a Debug folderben), akkor path is kell
fopen_s(&stream1, "A_proba_1.txt", "w");
fopen_s(&stream2, "A_proba_2.txt", "w");
fopen_s(&stream3, "A_proba_3.txt", "w");
fopen_s(&stream4, "A_proba_4.txt", "w");
fopen_s(&stream5, "A_proba_5.txt", "w");
fopen_s(&stream6, "A_proba_6.txt", "w");
//itt következik a kiindulási tartomány,
// a nyereghalmazt befoglaló tartomány:
long MMM = 10;
long NNN = 2;
long totalSize = static_cast<long>(pow((long double)(MMM + 1), 3) * pow((long double)(NNN + 1), 3));
clock_t startTime = clock();
#pragma omp parallel default(none) firstprivate(MMM, NNN, totalSize) shared(stream0, stream1, stream2, stream3, stream4, stream5, stream6)
{
#pragma omp for schedule(dynamic)
for (long n = 0; n < totalSize; n++) {
long i = n / (long)(pow((long double)(NNN + 1), 3) * pow((long double)(MMM + 1), 2)) % (MMM + 1);
long j = n / (long)(pow((long double)(NNN + 1), 3) * (MMM + 1)) % (MMM + 1);
long k = n / (long)pow((long double)(NNN + 1), 3) % (MMM + 1);
long ii = n / (long)pow((long double)(NNN + 1), 2) % (NNN + 1);
long jj = n / (NNN + 1) % (NNN + 1);
long kk = n % (NNN + 1);
/*for (long i = 0; i <= MMM; i++)
for (long j = 0; j <= MMM; j++)
for (long k = 0; k <= MMM; k++)
for (long ii = 0; ii <= NNN; ii++)
for (long jj = 0; jj <= NNN; jj++)
for (long kk = 0; kk <= NNN; kk++) {*/
double x[nmax + 1];
//az L5 pont környezetének kezdőfeltételei, ha az ábráról leolvasott értéket írok be, nem kell levonni a tkpx0-eket!!!
x[1] = x0 - 0.0002 + 0.0004* i / MMM - tkpx0; //ha az ábráról leolvasott értéket írok be, nem kell levonni a tkpx0-akat!!!
x[3] = yy0 - 0.0002 + 0.0004 * j / MMM - tkpy0;
x[5] = z0 - 0.0002 + 0.0004 * k / MMM - tkpz0;
x[2] = vx0 - 0.000015 + 0.00003* ii / NNN - tkpvx0;
x[4] = vy0 - 0.000015 + 0.00003* jj / NNN - tkpvy0;
x[6] = vz0 - 0.0000015 + 0.000003* kk / NNN - tkpvz0;
x[7] = nx0 - tkpx0; //Nap
x[8] = nvx0 - tkpvx0;
x[9] = nyy0 - tkpy0;
x[10] = nvy0 - tkpvy0;
x[11] = nz0 - tkpz0;
x[12] = nvy0 - tkpvz0;
x[13] = fx0 - tkpx0; //Föld
x[14] = fvx0 - tkpvx0;
x[15] = fyy0 - tkpy0;
x[16] = fvy0 - tkpvy0;
x[17] = fz0 - tkpz0;
x[18] = fvz0 - tkpvz0;
x[19] = hx0 - tkpx0; //Hold
x[20] = hvx0 - tkpvx0;
x[21] = hyy0 - tkpy0;
x[22] = hvy0 - tkpvy0;
x[23] = hz0 - tkpz0;
x[24] = hvz0 - tkpvz0;
Nyereghalmaz(/*hdc, rt,*/ i, j, k, ii, jj, kk, x, sizeof(x));
}
}
float elapsed = ((float)(clock() - startTime)) / CLOCKS_PER_SEC;
#ifndef _WIN32
elapsed = elapsed / (float)omp_get_num_procs();
#endif
printf("%.3f másodpercig futott.\n", elapsed);
fclose(stream0);
fclose(stream1);
fclose(stream2);
fclose(stream3);
fclose(stream4);
fclose(stream5);
fclose(stream6); // lezaras
}
int main(int argc, char* argv[]) {
#ifdef _WIN32
system("chcp 65001"); // Ez azért kell, hogy windowson is működjön az ékezetes betű kiírás
#endif
DrawContent();
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment