daily pastebin goal
76%
SHARE
TWEET

Untitled

a guest Jun 13th, 2018 54 in 16 hours
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. // geoph.cpp: определяет точку входа для консольного приложения.
  2. //
  3. #include "stdafx.h"
  4. #include <cmath>
  5.  
  6. const double f = 0.0000000000667; // TODO
  7.  
  8. double*** get_density(int grid_size) {
  9.     double*** density = new double**[grid_size];
  10.     for (int i = 0; i < grid_size; ++i) {
  11.         density[i] = new double*[grid_size];
  12.         for (int j = 0; j < grid_size; ++j) {
  13.             density[i][j] = new double[24];
  14.  
  15.             density[i][j][0] = 2000;
  16.  
  17.             double step = 30000.0 / grid_size;
  18.  
  19.             // 2 слой          
  20.             if (i < 14500 / step)
  21.                 density[i][j][1] = 2100;
  22.             else if (14500 / step <= i && i < 15500 / step)
  23.                 density[i][j][1] = 1900;
  24.             else
  25.                 density[i][j][1] = 2100;
  26.  
  27.             // 3 слой
  28.             if (i < 18000 / step)
  29.                 density[i][j][2] = 2200;
  30.             else if (18000 / step <= i && i < 20000 / step)
  31.                 density[i][j][2] = 2100;
  32.             else
  33.                 density[i][j][2] = 2200;
  34.  
  35.             // 4 слой
  36.             if (i < 1000 / step)
  37.                 density[i][j][3] = 2300;
  38.             else if (1000 / step <= i && i < 4000 / step)
  39.                 density[i][j][3] = 2250;
  40.             else if (4000 / step <= i && i < 26000 / step)
  41.                 density[i][j][3] = 2300;
  42.             else
  43.                 density[i][j][3] = 2250;
  44.  
  45.             for (int k = 4; k < 24; ++k) {
  46.                 density[i][j][k] = 2300 + (k - 4 + 1)*50;
  47.             }
  48.         }
  49.     }
  50.  
  51.     return density;
  52. }
  53.  
  54. double** get_initial_field_array(int grid_size) {
  55.     double** field = new double*[grid_size];
  56.     for (int i = 0; i < grid_size; ++i) {
  57.         field[i] = new double[grid_size];
  58.         for (int j = 0; j < grid_size; ++j) {
  59.             field[i][j] = 0;
  60.         }
  61.     }
  62.  
  63.     return field;
  64. }
  65.  
  66. void delete_density_array(double*** density_array, int grid_size) {
  67.     for (int i = 0; i < grid_size; ++i) {
  68.         for (int j = 0; j < 24; ++j) {
  69.             delete[] density_array[i][j];
  70.         }
  71.         delete[] density_array[i];
  72.     }
  73.  
  74.     delete[] density_array;
  75. }
  76.  
  77. void delete_field_array(double** field_array, double grid_size) {
  78.     for (int i = 0; i < grid_size; ++i) {
  79.         delete[] field_array[i];
  80.     }
  81.  
  82.     delete[] field_array;
  83. }
  84.  
  85. double get_norm(double x_coord, double y_coord, double z_coord) {
  86.     return sqrt(x_coord * x_coord + y_coord * y_coord + z_coord * z_coord);
  87. }
  88.  
  89. double get_v(double x_coord, double y_coord, double z_coord) {
  90.     auto norm = get_norm(x_coord, y_coord, z_coord);
  91.     return x_coord * log(y_coord + norm) + y_coord * log(x_coord + norm) - z_coord * atan(x_coord * y_coord / z_coord / norm);
  92. }
  93.  
  94. double get_density_or_zero(double*** density, int i, int j, int k, int grid_size, int layers_count)
  95. {
  96.     if (i == -1 || i == grid_size - 1)
  97.         return 0.0;
  98.     if (j == -1 || j == grid_size - 1)
  99.         return 0.0;
  100.     if (k == -1 || k == layers_count - 1)
  101.         return 0.0;
  102.  
  103.     return density[i][j][k];
  104. }
  105.  
  106. double get_reduced_density(double*** density, int i, int j, int k, int grid_size, int layers_count) {
  107.     double result_density = 0.0;
  108.  
  109.     result_density += get_density_or_zero(density, i - 1, j - 1, k - 1, grid_size, layers_count);
  110.     result_density -= get_density_or_zero(density, i - 1, j, k - 1, grid_size, layers_count); // density[i - 1][j][k - 1];
  111.     result_density -= get_density_or_zero(density, i, j - 1, k - 1, grid_size, layers_count); // density[i][j - 1][k - 1];
  112.     result_density += get_density_or_zero(density, i, j, k - 1, grid_size, layers_count); // density[i][j][k - 1];
  113.  
  114.     result_density -= get_density_or_zero(density, i - 1, j - 1, k, grid_size, layers_count); // density[i - 1][j - 1][k];
  115.     result_density += get_density_or_zero(density, i - 1, j, k, grid_size, layers_count); // density[i - 1][j][k];
  116.     result_density += get_density_or_zero(density, i, j - 1, k, grid_size, layers_count); // density[i][j - 1][k];
  117.     result_density -= get_density_or_zero(density, i, j, k, grid_size, layers_count); // density[i][j][k];
  118.  
  119.     return -result_density;
  120. }
  121.  
  122. double get_g(double*** density, int grid_size, int layers_count, double x_pos, double y_pos, double initial_shift, double step) {
  123.     double sum = 0.0;
  124.     double x = initial_shift;
  125.     double y;
  126.     for (int i = 0; i < grid_size + 1; ++i) {
  127.         y = initial_shift;
  128.         for (int j = 0; j < grid_size + 1; ++j) {
  129.             for (int k = 0; k < layers_count; ++k) {
  130.                 sum += get_reduced_density(density, i, j, k, grid_size, layers_count) * get_v(x - x_pos, y - y_pos, -(k + 1)*1000);// TODO: надо определить шаг по z
  131.                 //sum += get_reduced_density(density, i, j, k, grid_size, layers_count) * get_v(x_pos - x, y_pos - y, (k + 1) * 1000);
  132.                 //sum += get_reduced_density(density, i, j, k, grid_size, layers_count) * get_v(x_pos - x, y_pos - y, -(k + 1)*1000);
  133.             }
  134.             y += step;
  135.         }
  136.         x += step;
  137.     }
  138.  
  139.     return f * sum;
  140. }
  141.  
  142. int main() {
  143.     int grid_size = 60;
  144.     auto density = get_density(grid_size + 1);
  145.     double step = 30000.0 / grid_size;
  146.     int pl_size = 2 * grid_size;
  147.     auto field = get_initial_field_array(pl_size + 1);
  148.     double initial_shift = grid_size / 2.0 * step;
  149.     FILE *file;
  150.     file = fopen("file.dat", "w");
  151.     if (file != NULL) {
  152.         for (int i = 0; i < pl_size; ++i) {
  153.             for (int j = 0; j < pl_size; ++j) {
  154.                 field[i][j] = get_g(density, grid_size, 24, i*step, j*step, initial_shift, step); // TODO: layers_count
  155.                 fprintf(file, "%lf %lf %lf\r\n", i*step, j*step, field[i][j]);
  156.             }
  157.         }
  158.     }
  159.  
  160.     delete_density_array(density, grid_size);
  161.     delete_field_array(field, pl_size);
  162.     // TODO: вывод массива
  163.     return 0;
  164. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top