780 lines
24 KiB
C++
780 lines
24 KiB
C++
/********************************************************
|
|
* ██████╗ ██████╗████████╗██╗
|
|
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
|
* ██║ ███╗██║ ██║ ██║
|
|
* ██║ ██║██║ ██║ ██║
|
|
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
|
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
|
* Geophysical Computational Tools & Library (GCTL)
|
|
*
|
|
* Copyright (c) 2022 Yi Zhang (yizhang-geo@zju.edu.cn)
|
|
*
|
|
* GCTL is distributed under a dual licensing scheme. You can redistribute
|
|
* it and/or modify it under the terms of the GNU Lesser General Public
|
|
* License as published by the Free Software Foundation, either version 2
|
|
* of the License, or (at your option) any later version. You should have
|
|
* received a copy of the GNU Lesser General Public License along with this
|
|
* program. If not, see <http://www.gnu.org/licenses/>.
|
|
*
|
|
* If the terms and conditions of the LGPL v.2. would prevent you from using
|
|
* the GCTL, please consider the option to obtain a commercial license for a
|
|
* fee. These licenses are offered by the GCTL's original author. As a rule,
|
|
* licenses are provided "as-is", unlimited in time for a one time fee. Please
|
|
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
|
|
* to include some description of your company and the realm of its activities.
|
|
* Also add information on how to contact you by electronic and paper mail.
|
|
******************************************************/
|
|
|
|
#ifndef _FUNC_H
|
|
#define _FUNC_H
|
|
|
|
// add gctl head files
|
|
#include "gctl/core.h"
|
|
#include "gctl/utility.h"
|
|
#include "gctl/math/legendre.h"
|
|
#include "gctl/math/interpolate.h"
|
|
#include "gctl/io.h"
|
|
|
|
#define NORMAL_GRAVITY 9.80665 // m/s^2
|
|
|
|
enum target_type_e
|
|
{
|
|
Null,
|
|
Topography,
|
|
GravDisturbance,
|
|
GravDisturbanceLon,
|
|
GravDisturbanceLat,
|
|
GravAnomaly,
|
|
GravPotential,
|
|
HeightAnomaly,
|
|
RadialGravityGradient,
|
|
};
|
|
|
|
struct spoint : public gctl::point3ds
|
|
{
|
|
double ref, alti, val;
|
|
spoint()
|
|
{
|
|
lon = lat = ref = alti = rad = val = GCTL_BDL_MAX;
|
|
}
|
|
void info()
|
|
{
|
|
std::cout << std::setprecision(16) << lon << " " << lat << " " << ref << " " << alti << " " << val << std::endl;
|
|
}
|
|
};
|
|
typedef std::vector<spoint> sphArray;
|
|
|
|
// 命令规则 n为阶数 m为次数
|
|
class shc2xyz
|
|
{
|
|
public:
|
|
shc2xyz(){}
|
|
~shc2xyz(){}
|
|
int readSHC(const char*,const char*,const char*, std::string jp); //读入球谐系数
|
|
int initTargetType(const char*); //设置计算类型
|
|
int initMatrix(const char*,const char*,const char*,const char*); //初始化相关的矩阵大小
|
|
int initObs(const char*,const char*,const char*); //初始化观测点 如只有范围参数则只初始化经纬度位置
|
|
int relocateAltitude(const char*); //根据输入文件重新确定计算高程
|
|
int outObs(const char*); //输出计算结果 如果有文件指定的位置则插值
|
|
int calSolution(); //计算球谐结果 同一高程观测值
|
|
int calSolution2(const char*); //计算不同高程的观测值
|
|
private:
|
|
gctl::array<gctl::array<double>> Anm;
|
|
gctl::array<gctl::array<double>> Bnm;
|
|
gctl::array<gctl::array<double>> Pnm; //伴随勒让德函数系数 这个函数只和观测位置的纬度/余纬度相关 同一纬度只需要计算一次
|
|
gctl::matrix<double> mCos; //不同次数cos函数值 这个值只和观测位置的经度相关 行数为不同经度位置 列数为不同次数 矩阵维度即为经度个数*阶次 一般估算在1000*1000级别
|
|
gctl::matrix<double> mSin; //不同次数sin函数值 其他与上同
|
|
gctl::array<gctl::array<double>> coff_S; //球谐系数sin参数
|
|
gctl::array<gctl::array<double>> coff_C; //球谐系数cos参数
|
|
gctl::array<gctl::array<double>> multi_array; //乘子矩阵
|
|
sphArray obsPoint; //计算地形是的观测位置 即计算半径值
|
|
sphArray outPoint; //输出计算值
|
|
|
|
gctl::legendre_norm_e norSum;
|
|
double GM,R; //球谐系数中重力常数与质量的乘积 单位为SI标准 g 与 m
|
|
double multi_factor; // 乘子系数
|
|
int NN_size; //系数矩阵大小
|
|
int lon_size,lat_size;
|
|
double refr,refR,altitude;
|
|
double lonmin,lonmax,dlon;
|
|
double latmin,latmax,dlat;
|
|
|
|
target_type_e tar_type_;
|
|
};
|
|
|
|
//读取球谐参数文件 文件名 起止阶次 列序列
|
|
int shc2xyz::readSHC(const char* filename, const char* para, const char* orders, std::string jp)
|
|
{
|
|
std::ifstream infile;
|
|
gctl::open_infile(infile,filename, "");
|
|
|
|
int n_start,m_start,n_end,m_end;
|
|
if (4 != sscanf(para,"%d/%d/%d/%d",&n_start,&m_start,&n_end,&m_end))
|
|
{
|
|
std::cerr << GCTL_BOLDRED << "error ==> " << GCTL_RESET << "wrong parameter: " << para << std::endl;
|
|
return -1;
|
|
}
|
|
//识别列次序
|
|
int order[4];
|
|
if (4 != sscanf(orders,"%d,%d,%d,%d",&order[0],&order[1],&order[2],&order[3]))
|
|
{
|
|
std::cerr << GCTL_BOLDRED << "error ==> " << GCTL_RESET << "wrong parameter: " << orders << std::endl;
|
|
return -1;
|
|
}
|
|
|
|
//按照最大阶数初始化下半三角矩阵
|
|
NN_size = n_end + 1;
|
|
//对于二维vector来说 对行初始化的时候需要使用resize 而对于列的初始化而言使用reserve效率更高
|
|
coff_C.resize(NN_size);
|
|
coff_S.resize(NN_size);
|
|
for (int i = 0; i < NN_size; i++)
|
|
{
|
|
coff_C[i].resize(i+1,0.0);
|
|
coff_S[i].resize(i+1,0.0);
|
|
}
|
|
|
|
int jp_num = atoi(jp.c_str());
|
|
int n,m; //行列号
|
|
std::string temp_d;
|
|
double temp_c,temp_s;
|
|
std::vector<std::string> temp_row; temp_row.reserve(100); //出现初始化100个double的空间 这样读文件更快
|
|
std::string temp_str;
|
|
std::stringstream temp_ss;
|
|
while (getline(infile, temp_str))
|
|
{
|
|
if (jp_num > 0) {jp_num--; continue;}
|
|
if (temp_str[0] == '#') continue;
|
|
|
|
//gctl::parse_string_to_vector(temp_str, ' ', temp_row);
|
|
|
|
if (*(temp_str.begin()) == '#') continue;
|
|
if (!temp_row.empty()) temp_row.clear();
|
|
temp_ss.str("");
|
|
temp_ss.clear();
|
|
temp_ss << temp_str;
|
|
while (temp_ss >> temp_d)
|
|
{
|
|
temp_row.push_back(temp_d);
|
|
}
|
|
|
|
n = atoi(temp_row[order[0]].c_str());
|
|
m = atoi(temp_row[order[1]].c_str());
|
|
temp_c = gctl::str2double(temp_row[order[2]]);
|
|
temp_s = gctl::str2double(temp_row[order[3]]);
|
|
|
|
if (n >= n_start && n <= n_end && m >= m_start && m <= m_end)
|
|
{
|
|
coff_C[n][m] = temp_c;
|
|
coff_S[n][m] = temp_s;
|
|
}
|
|
}
|
|
infile.close();
|
|
return 0;
|
|
}
|
|
|
|
int shc2xyz::initObs(const char* r_para,const char* i_para,const char* refsys)
|
|
{
|
|
//解析参考球
|
|
if (!strcmp(refsys,"NULL"))
|
|
{
|
|
refr = refR = 0.0;
|
|
}
|
|
else if (!strcmp(refsys,"WGS84"))
|
|
{
|
|
refr = GCTL_WGS84_PoleRadius;
|
|
refR = GCTL_WGS84_EquatorRadius;
|
|
}
|
|
else if (!strcmp(refsys,"Earth"))
|
|
{
|
|
refr = GCTL_Earth_Radius;
|
|
refR = GCTL_Earth_Radius;
|
|
}
|
|
else if (!strcmp(refsys,"Moon"))
|
|
{
|
|
refr = GCTL_Moon_Radius;
|
|
refR = GCTL_Moon_Radius;
|
|
}
|
|
else if (!strcmp(refsys,"Mars"))
|
|
{
|
|
refr = GCTL_Mars_Radius;
|
|
refR = GCTL_Mars_Radius;
|
|
}
|
|
else if (2 != sscanf(refsys,"%lf/%lf",&refr,&refR))
|
|
{
|
|
std::cerr << GCTL_BOLDRED << "error ==> " << GCTL_RESET << "wrong parameter: " << refsys << std::endl;
|
|
return -1;
|
|
}
|
|
|
|
//解析经纬度范围 按规则网络初始化观测点位置
|
|
if (5 != sscanf(r_para,"%lf/%lf/%lf/%lf/%lf",&lonmin,&lonmax,&latmin,&latmax,&altitude))
|
|
{
|
|
if (4 != sscanf(r_para,"%lf/%lf/%lf/%lf",&lonmin,&lonmax,&latmin,&latmax))
|
|
{
|
|
std::cerr << GCTL_BOLDRED << "error ==> " << GCTL_RESET << "wrong parameter: " << r_para << std::endl;
|
|
return -1;
|
|
}
|
|
else altitude = 0.0;
|
|
}
|
|
|
|
//解析间隔
|
|
if (2 != sscanf(i_para,"%lf/%lf",&dlon,&dlat))
|
|
{
|
|
std::cerr << GCTL_BOLDRED << "error ==> " << GCTL_RESET << "wrong parameter: " << i_para << std::endl;
|
|
return -1;
|
|
}
|
|
|
|
spoint temp_spoint;
|
|
double lon,lat;
|
|
lon_size = round((lonmax-lonmin)/dlon) + 1;
|
|
lat_size = round((latmax-latmin)/dlat) + 1;
|
|
obsPoint.reserve(lon_size*lat_size);
|
|
for (int i = 0; i < lat_size; i++)
|
|
{
|
|
for (int j = 0; j < lon_size; j++)
|
|
{
|
|
lat = latmin + i*dlat; lon = lonmin + j*dlon;
|
|
temp_spoint.lon = lon; temp_spoint.lat = lat;
|
|
temp_spoint.ref = gctl::ellipse_radius_2d(refR, refr, temp_spoint.lat*GCTL_Pi/180.0);
|
|
temp_spoint.alti = altitude;
|
|
temp_spoint.rad = temp_spoint.ref + temp_spoint.alti;
|
|
temp_spoint.val = 0.0;
|
|
obsPoint.push_back(temp_spoint);
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
int shc2xyz::relocateAltitude(const char* filepara)
|
|
{
|
|
char filename[1024];
|
|
int orders[3] = {0,1,2}; //默认的读入的数据列为前三列
|
|
|
|
if(!strcmp(filepara,"NULL")) return 0;
|
|
|
|
//解析文件名中是否含有+d标示 如果有则将+d以前解释为filename 之后为需要读入的数据列 默认为逗号分隔
|
|
//否则将filepara赋值为filename
|
|
if (4 != sscanf(filepara,"%[^+]+d%d,%d,%d",filename,&orders[0],&orders[1],&orders[2]))
|
|
strcpy(filename,filepara);
|
|
|
|
std::ifstream infile;
|
|
gctl::open_infile(infile,filename,"");
|
|
|
|
int numM,numN,tempM,tempN;
|
|
std::string temp_str;
|
|
std::stringstream temp_ss;
|
|
double temp_d,temp_lon,temp_lat,temp_alti;
|
|
std::vector<double> temp_row;
|
|
|
|
numM = floor((latmax-latmin)/dlat)+1;
|
|
numN = floor((lonmax-lonmin)/dlon)+1;
|
|
while(getline(infile,temp_str))
|
|
{
|
|
if (*(temp_str.begin()) == '#') continue;
|
|
|
|
temp_ss.str(""); temp_ss.clear(); temp_ss << temp_str;
|
|
if(!temp_row.empty()) temp_row.clear();
|
|
while(temp_ss >> temp_d)
|
|
temp_row.push_back(temp_d);
|
|
|
|
temp_lon = temp_row[orders[0]];
|
|
temp_lat = temp_row[orders[1]];
|
|
temp_alti = temp_row[orders[2]];
|
|
|
|
tempM = round((temp_lat-latmin)/dlat);
|
|
tempN = round((temp_lon-lonmin)/dlon);
|
|
obsPoint[tempM*numN+tempN].alti = temp_alti;
|
|
obsPoint[tempM*numN+tempN].rad = obsPoint[tempM*numN+tempN].ref + temp_alti;
|
|
}
|
|
infile.close();
|
|
return 0;
|
|
}
|
|
|
|
int shc2xyz::initTargetType(const char* type)
|
|
{
|
|
if(!strcmp(type,"n")) //topography
|
|
tar_type_ = Null;
|
|
else if(!strcmp(type,"t")) //topography
|
|
tar_type_ = Topography;
|
|
else if (!strcmp(type,"d") || !strcmp(type,"g"))
|
|
tar_type_ = GravAnomaly;
|
|
else if (!strcmp(type,"r"))
|
|
tar_type_ = RadialGravityGradient;
|
|
else if (!strcmp(type,"p"))
|
|
tar_type_ = GravPotential;
|
|
else if (!strcmp(type,"h"))
|
|
tar_type_ = HeightAnomaly;
|
|
else if (!strcmp(type,"dlon"))
|
|
tar_type_ = GravDisturbanceLon;
|
|
else if (!strcmp(type,"dlat"))
|
|
tar_type_ = GravDisturbanceLat;
|
|
else
|
|
{
|
|
std::cerr << GCTL_BOLDRED << "error ==> " << GCTL_RESET << "unknown calculation type of " << type << std::endl;
|
|
return -1;
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
//初始化矩阵
|
|
int shc2xyz::initMatrix(const char* type,const char* para,const char* norType,const char* zfile)
|
|
{
|
|
//初始化GM与R
|
|
if (strcmp(para,"NULL")) //如果para不为NULL则识别参数 否则将GM与R初始化为MAX_BDL
|
|
{
|
|
if (2 != sscanf(para,"%lf/%lf",&GM,&R))
|
|
{
|
|
std::cerr << GCTL_BOLDRED << "error ==> " << GCTL_RESET << "wrong parameter: " << para << std::endl;
|
|
return -1;
|
|
}
|
|
}
|
|
else GM = R = GCTL_BDL_MAX;
|
|
//初始化归一化类型
|
|
if (!strcmp(norType,"g")) norSum = gctl::Pi4;
|
|
else if (!strcmp(norType,"m")) norSum = gctl::One;
|
|
else
|
|
{
|
|
std::cerr << GCTL_BOLDRED << "error ==> " << GCTL_RESET << "wrong parameter: " << norType << std::endl;
|
|
return -1;
|
|
}
|
|
//初始化伴随勒让德函数矩阵
|
|
Pnm.resize(NN_size);
|
|
for (int i = 0; i < NN_size; i++)
|
|
Pnm.at(i).resize(i+1,0.0);
|
|
//初始化sin和cos矩阵
|
|
mCos.resize(lon_size, NN_size);
|
|
mSin.resize(lon_size, NN_size);
|
|
//计算mCos和mSin的值
|
|
int i,j;
|
|
double lon;
|
|
#pragma omp parallel for private(i,j,lon) schedule(guided)
|
|
for (i = 0; i < lon_size; i++)
|
|
{
|
|
lon = lonmin + i*dlon;
|
|
for (j = 0; j < NN_size; j++)
|
|
{
|
|
mCos[i][j] = cos(j*lon*GCTL_Pi/180.0);
|
|
mSin[i][j] = sin(j*lon*GCTL_Pi/180.0);
|
|
}
|
|
}
|
|
//计算勒让德函数系数
|
|
gctl::get_a_nm_array(NN_size, Anm);
|
|
gctl::get_b_nm_array(NN_size, Bnm);
|
|
//计算乘子参数
|
|
if(!strcmp(type,"n")) // null
|
|
multi_factor = 1.0;
|
|
else if(!strcmp(type,"t")) //topography
|
|
multi_factor = 1.0;
|
|
else if (!strcmp(type,"d") || !strcmp(type,"g")) //gravity disturbance
|
|
multi_factor = 1e+5*GM/(R*R);
|
|
else if (!strcmp(type,"r")) //gravity disturbance
|
|
multi_factor = 1e+9*GM/(R*R);
|
|
else if (!strcmp(type,"p"))
|
|
multi_factor = 1e+5*GM/R;
|
|
else if (!strcmp(type,"h"))
|
|
multi_factor = 1e+5*GM/(R*NORMAL_GRAVITY*1e+5);
|
|
else if (!strcmp(type,"dlat") || !strcmp(type,"dlon"))
|
|
multi_factor = 1e+5*GM/(R*R);
|
|
else
|
|
{
|
|
std::cerr << GCTL_BOLDRED << "error ==> " << GCTL_RESET << "unknown calculation type of " << type << std::endl;
|
|
return -1;
|
|
}
|
|
//初始化乘子矩阵
|
|
multi_array.resize(lat_size);
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
multi_array[i].resize(NN_size,1.0); //初始化乘子矩阵为1 适用于地形等直接计算的类型
|
|
}
|
|
//如果计算高程不在同一高程 则不能使用multi_array 同时应该使用calSolution2()函数
|
|
if (strcmp(zfile,"NULL"))
|
|
return 0;
|
|
//如果计算类型不是地形等直接计算类型则需要检验-g选项是否已经设置
|
|
if (strcmp(type,"t"))
|
|
{
|
|
if (GM == GCTL_BDL_MAX || R == GCTL_BDL_MAX)
|
|
{
|
|
std::cerr << GCTL_BOLDRED << "error ==> " << GCTL_RESET << "-g option must be set for gravitational calculation" << std::endl;
|
|
return -1;
|
|
}
|
|
}
|
|
//根据不同类型计算乘子参数和乘子矩阵
|
|
if (!strcmp(type,"d")) //gravity disturbance
|
|
{
|
|
//#pragma omp parallel for private(i,j) shared(R,lon_size) schedule(guided)
|
|
#pragma omp parallel for private(i,j) schedule(guided)
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
for (j = 0; j < NN_size; j++)
|
|
{
|
|
multi_array[i][j] = pow(R/obsPoint[i*lon_size].rad,j+2)*(j+1);
|
|
}
|
|
}
|
|
}
|
|
else if (!strcmp(type,"r")) //gravity gradient
|
|
{
|
|
//#pragma omp parallel for private(i,j) shared(R,lon_size) schedule(guided)
|
|
#pragma omp parallel for private(i,j) schedule(guided)
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
for (j = 0; j < NN_size; j++)
|
|
{
|
|
multi_array[i][j] = pow(R/obsPoint[i*lon_size].rad,j+2)*(j+1)*(j+2)/obsPoint[i*lon_size].rad;
|
|
}
|
|
}
|
|
}
|
|
else if (!strcmp(type,"g")) //gravity anomaly
|
|
{
|
|
//#pragma omp parallel for private(i,j) shared(R,lon_size) schedule(guided)
|
|
#pragma omp parallel for private(i,j) schedule(guided)
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
for (j = 0; j < NN_size; j++)
|
|
{
|
|
multi_array[i][j] = pow(R/obsPoint[i*lon_size].rad,j+2)*(j-1);
|
|
}
|
|
}
|
|
}
|
|
else if (!strcmp(type,"p")) //geo-potential
|
|
{
|
|
//#pragma omp parallel for private(i,j) shared(R,lon_size) schedule(guided)
|
|
#pragma omp parallel for private(i,j) schedule(guided)
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
for (j = 0; j < NN_size; j++)
|
|
{
|
|
multi_array[i][j] = pow(R/obsPoint[i*lon_size].rad,j+1);
|
|
}
|
|
}
|
|
}
|
|
else if (!strcmp(type,"dlon") || !strcmp(type,"dlat"))
|
|
{
|
|
//#pragma omp parallel for private(i,j) shared(R,lon_size) schedule(guided)
|
|
#pragma omp parallel for private(i,j) schedule(guided)
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
for (j = 0; j < NN_size; j++)
|
|
{
|
|
multi_array[i][j] = pow(R/obsPoint[i*lon_size].rad,j+2);
|
|
}
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
int shc2xyz::outObs(const char* filename)
|
|
{
|
|
if (!strcmp(filename,"NULL")) //没有输入文件 直接输出规则网计算结果
|
|
{
|
|
if (tar_type_ == Null)
|
|
{
|
|
std::cout << "# NaN value = 1e+30" << std::endl;
|
|
std::cout << "# lon(deg) lat(deg) data(unknown)" << std::endl;
|
|
for (int i = 0; i < obsPoint.size(); i++)
|
|
{
|
|
std::cout << std::setprecision(16) << obsPoint[i].lon << " " << obsPoint[i].lat << " " << obsPoint[i].val << std::endl;
|
|
}
|
|
}
|
|
else if (tar_type_ == Topography)
|
|
{
|
|
std::cout << "# NaN value = 1e+30" << std::endl;
|
|
std::cout << "# topography = radius - reference" << std::endl;
|
|
std::cout << "# lon(deg) lat(deg) reference-radius(m) topography(m)" << std::endl;
|
|
for (int i = 0; i < obsPoint.size(); i++)
|
|
{
|
|
std::cout << std::setprecision(16) << obsPoint[i].lon << " " << obsPoint[i].lat << " "
|
|
<< obsPoint[i].ref << " " << obsPoint[i].val - obsPoint[i].ref << std::endl;
|
|
}
|
|
}
|
|
else if (tar_type_ == HeightAnomaly)
|
|
{
|
|
std::cout << "# NaN value = 1e+30" << std::endl;
|
|
std::cout << "# Normal Gravity = " << NORMAL_GRAVITY << " m/s^2" << std::endl;
|
|
std::cout << "# lon(deg) lat(deg) reference-radius(m) height-anomaly(m)" << std::endl;
|
|
for (int i = 0; i < obsPoint.size(); i++)
|
|
{
|
|
obsPoint[i].info();
|
|
}
|
|
}
|
|
else
|
|
{
|
|
std::cout << "# NaN value = 1e+30" << std::endl;
|
|
std::cout << "# lon(deg) lat(deg) reference-radius(m) altitude(m) grav-field(mGal|Eo)" << std::endl;
|
|
for (int i = 0; i < obsPoint.size(); i++)
|
|
{
|
|
obsPoint[i].info();
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
std::ifstream infile;
|
|
gctl::open_infile(infile,filename,"");
|
|
|
|
spoint temp_sp;
|
|
std::string temp_str;
|
|
std::stringstream temp_ss;
|
|
while (getline(infile,temp_str))
|
|
{
|
|
if(*(temp_str.begin()) == '#') continue;
|
|
|
|
temp_ss.str("");
|
|
temp_ss.clear();
|
|
temp_ss << temp_str;
|
|
temp_ss >> temp_sp.lon >> temp_sp.lat;
|
|
temp_sp.ref = gctl::ellipse_radius_2d(refR, refr, temp_sp.lat*GCTL_Pi/180.0);
|
|
temp_sp.alti = altitude;
|
|
outPoint.push_back(temp_sp);
|
|
}
|
|
infile.close();
|
|
|
|
int numM,numN,tempM,tempN;
|
|
double lon1,lon2,lat1,lat2;
|
|
|
|
numM = floor((latmax-latmin)/dlat)+1;
|
|
numN = floor((lonmax-lonmin)/dlon)+1;
|
|
for (int i = 0; i < outPoint.size(); i++)
|
|
{
|
|
tempM = floor((outPoint[i].lat-latmin)/dlat);
|
|
tempN = floor((outPoint[i].lon-lonmin)/dlon);
|
|
|
|
if (tempM == (numM-1))
|
|
tempM -= 1;
|
|
if (tempN == (numN-1))
|
|
tempN -= 1;
|
|
|
|
if (tempM >= 0 && tempN >= 0 && tempM <= numM-2 && tempN <= numN-2)
|
|
{
|
|
lon1 = lonmin+tempN*dlon;
|
|
lon2 = lonmin+(tempN+1)*dlon;
|
|
lat1 = latmin+tempM*dlat;
|
|
lat2 = latmin+(tempM+1)*dlat;
|
|
outPoint[i].val = gctl::sph_linear_interpolate_deg(lat1,lat2,lon1,lon2,
|
|
outPoint[i].lat,outPoint[i].lon,
|
|
obsPoint[tempM*numN+tempN].val,
|
|
obsPoint[tempM*numN+tempN+1].val,
|
|
obsPoint[(tempM+1)*numN+tempN].val,
|
|
obsPoint[(tempM+1)*numN+tempN+1].val);
|
|
}
|
|
}
|
|
|
|
if (tar_type_ == Null)
|
|
{
|
|
std::cout << "# NaN value = 1e+30" << std::endl;
|
|
std::cout << "# lon(deg) lat(deg) data(unknown)" << std::endl;
|
|
for (int i = 0; i < outPoint.size(); i++)
|
|
{
|
|
std::cout << std::setprecision(16) << outPoint[i].lon << " " << outPoint[i].lat << " " << outPoint[i].val << std::endl;
|
|
}
|
|
}
|
|
else if (tar_type_ == Topography)
|
|
{
|
|
std::cout << "# NaN value = 1e+30" << std::endl;
|
|
std::cout << "# topography = radius - reference" << std::endl;
|
|
std::cout << "# lon(deg) lat(deg) reference-radius(m) topography(m)" << std::endl;
|
|
for (int i = 0; i < outPoint.size(); i++)
|
|
{
|
|
std::cout << std::setprecision(16) << outPoint[i].lon << " " << outPoint[i].lat << " "
|
|
<< outPoint[i].ref << " " << outPoint[i].val - outPoint[i].ref << std::endl;
|
|
}
|
|
}
|
|
else if (tar_type_ == HeightAnomaly)
|
|
{
|
|
std::cout << "# NaN value = 1e+30" << std::endl;
|
|
std::cout << "# Normal Gravity = " << NORMAL_GRAVITY << " m/s^2" << std::endl;
|
|
std::cout << "# lon(deg) lat(deg) reference-radius(m) height-anomaly(m)" << std::endl;
|
|
for (int i = 0; i < outPoint.size(); i++)
|
|
{
|
|
outPoint[i].info();
|
|
}
|
|
}
|
|
else
|
|
{
|
|
std::cout << "# NaN value = 1e+30" << std::endl;
|
|
std::cout << "# lon(deg) lat(deg) reference-radius(m) altitude(m) grav-field(mGal|Eo)" << std::endl;
|
|
for (int i = 0; i < outPoint.size(); i++)
|
|
{
|
|
outPoint[i].info();
|
|
}
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
int shc2xyz::calSolution()
|
|
{
|
|
//计算
|
|
int i,j,n,m;
|
|
double temp_d,lat;
|
|
|
|
gctl::progress_bar bar(lat_size,"Process");
|
|
|
|
if (tar_type_ == GravDisturbanceLon)
|
|
{
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
bar.progressed(i);
|
|
lat = latmin + dlat*i;
|
|
//计算伴随勒让德函数 对于同一个纬度只需要计算一次
|
|
gctl::nalf_sfcm(Pnm,Anm,Bnm,NN_size,90.0-lat,norSum);
|
|
//这里可以使用并行加速计算外层循环 内层计算因为是递归计算因此不能并行
|
|
//一种并行方案更快一些
|
|
#pragma omp parallel for private(j,n,m,temp_d) shared(i) schedule(guided)
|
|
for (j = 0; j < lon_size; j++)
|
|
{
|
|
temp_d = 0;
|
|
for (n = 0; n < NN_size; n++)
|
|
{
|
|
for (m = 0; m < n+1; m++)
|
|
{
|
|
temp_d += (1.0/cos(lat*GCTL_Pi/180.0))*multi_array[i][n]*m*Pnm[n][m]*(coff_S[n][m]*mCos[j][m] - coff_C[n][m]*mSin[j][m]);
|
|
}
|
|
}
|
|
obsPoint[i*lon_size+j].val = multi_factor*temp_d;
|
|
}
|
|
}
|
|
}
|
|
else if (tar_type_ == GravDisturbanceLat)
|
|
{
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
bar.progressed(i);
|
|
lat = latmin + dlat*i;
|
|
//计算伴随勒让德函数 对于同一个纬度只需要计算一次
|
|
gctl::nalf_sfcm(Pnm,Anm,Bnm,NN_size,90.0-lat,norSum,true);
|
|
//这里可以使用并行加速计算外层循环 内层计算因为是递归计算因此不能并行
|
|
//一种并行方案更快一些
|
|
#pragma omp parallel for private(j,n,m,temp_d) shared(i) schedule(guided)
|
|
for (j = 0; j < lon_size; j++)
|
|
{
|
|
temp_d = 0;
|
|
for (n = 0; n < NN_size; n++)
|
|
{
|
|
for (m = 0; m < n+1; m++)
|
|
{
|
|
temp_d += multi_array[i][n]*Pnm[n][m]*(coff_C[n][m]*mCos[j][m]+coff_S[n][m]*mSin[j][m]);
|
|
}
|
|
}
|
|
obsPoint[i*lon_size+j].val = multi_factor*temp_d;
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
bar.progressed(i);
|
|
lat = latmin + dlat*i;
|
|
//计算伴随勒让德函数 对于同一个纬度只需要计算一次
|
|
gctl::nalf_sfcm(Pnm,Anm,Bnm,NN_size,90.0-lat,norSum,false);
|
|
//这里可以使用并行加速计算外层循环 内层计算因为是递归计算因此不能并行
|
|
//一种并行方案更快一些
|
|
#pragma omp parallel for private(j,n,m,temp_d) shared(i) schedule(guided)
|
|
for (j = 0; j < lon_size; j++)
|
|
{
|
|
temp_d = 0;
|
|
for (n = 0; n < NN_size; n++)
|
|
{
|
|
for (m = 0; m < n+1; m++)
|
|
{
|
|
temp_d += multi_array[i][n]*Pnm[n][m]*(coff_C[n][m]*mCos[j][m]+coff_S[n][m]*mSin[j][m]);
|
|
}
|
|
}
|
|
obsPoint[i*lon_size+j].val = multi_factor*temp_d;
|
|
}
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
int shc2xyz::calSolution2(const char* type)
|
|
{
|
|
//计算
|
|
int i,j,n,m;
|
|
double temp_d,lat;
|
|
|
|
gctl::progress_bar bar(lat_size,"Process");
|
|
if (!strcmp(type,"d"))
|
|
{
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
bar.progressed(i);
|
|
lat = latmin + dlat*i;
|
|
//计算伴随勒让德函数 对于同一个纬度只需要计算一次
|
|
gctl::nalf_sfcm(Pnm,Anm,Bnm,NN_size,90.0-lat,norSum);
|
|
//这里可以使用并行加速计算外层循环 内层计算因为是递归计算因此不能并行
|
|
//一种并行方案更快一些
|
|
//#pragma omp parallel for private(j,n,m,temp_d) shared(i,multi_factor,lon_size) schedule(guided)
|
|
#pragma omp parallel for private(j,n,m,temp_d) shared(i) schedule(guided)
|
|
for (j = 0; j < lon_size; j++)
|
|
{
|
|
temp_d = 0;
|
|
for (n = 0; n < NN_size; n++)
|
|
{
|
|
for (m = 0; m < n+1; m++)
|
|
{
|
|
//pow(R/obsPoint[i*lon_size+j].rad,n+2)*(n+1)
|
|
temp_d += pow(R/obsPoint[i*lon_size+j].rad,n+2)*(n+1)*Pnm[n][m]*(coff_C[n][m]*mCos[j][m]+coff_S[n][m]*mSin[j][m]);
|
|
}
|
|
}
|
|
obsPoint[i*lon_size+j].val = multi_factor*temp_d;
|
|
}
|
|
}
|
|
}
|
|
else if (!strcmp(type,"g"))
|
|
{
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
bar.progressed(i);
|
|
lat = latmin + dlat*i;
|
|
//计算伴随勒让德函数 对于同一个纬度只需要计算一次
|
|
gctl::nalf_sfcm(Pnm,Anm,Bnm,NN_size,90.0-lat,norSum);
|
|
//这里可以使用并行加速计算外层循环 内层计算因为是递归计算因此不能并行
|
|
//一种并行方案更快一些
|
|
//#pragma omp parallel for private(j,n,m,temp_d) shared(i,multi_factor,lon_size) schedule(guided)
|
|
#pragma omp parallel for private(j,n,m,temp_d) shared(i) schedule(guided)
|
|
for (j = 0; j < lon_size; j++)
|
|
{
|
|
temp_d = 0;
|
|
for (n = 0; n < NN_size; n++)
|
|
{
|
|
for (m = 0; m < n+1; m++)
|
|
{
|
|
//pow(R/obsPoint[i*lon_size+j].rad,n+2)*(n-1)
|
|
temp_d += pow(R/obsPoint[i*lon_size+j].rad,n+2)*(n-1)*Pnm[n][m]*(coff_C[n][m]*mCos[j][m]+coff_S[n][m]*mSin[j][m]);
|
|
}
|
|
}
|
|
obsPoint[i*lon_size+j].val = multi_factor*temp_d;
|
|
}
|
|
}
|
|
}
|
|
else if (!strcmp(type,"p"))
|
|
{
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
bar.progressed(i);
|
|
lat = latmin + dlat*i;
|
|
//计算伴随勒让德函数 对于同一个纬度只需要计算一次
|
|
gctl::nalf_sfcm(Pnm,Anm,Bnm,NN_size,90.0-lat,norSum);
|
|
//这里可以使用并行加速计算外层循环 内层计算因为是递归计算因此不能并行
|
|
//一种并行方案更快一些
|
|
//#pragma omp parallel for private(j,n,m,temp_d) shared(i,multi_factor,lon_size) schedule(guided)
|
|
#pragma omp parallel for private(j,n,m,temp_d) shared(i) schedule(guided)
|
|
for (j = 0; j < lon_size; j++)
|
|
{
|
|
temp_d = 0;
|
|
for (n = 0; n < NN_size; n++)
|
|
{
|
|
for (m = 0; m < n+1; m++)
|
|
{
|
|
//pow(R/obsPoint[i*lon_size+j].rad,n+1)
|
|
temp_d += pow(R/obsPoint[i*lon_size+j].rad,n+1)*Pnm[n][m]*(coff_C[n][m]*mCos[j][m]+coff_S[n][m]*mSin[j][m]);
|
|
}
|
|
}
|
|
obsPoint[i*lon_size+j].val = multi_factor*temp_d;
|
|
}
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
#endif //_FUNC_H
|