432 lines
11 KiB
C++
432 lines
11 KiB
C++
#ifndef _FUNC_H
|
|
#define _FUNC_H
|
|
#include "sysDefine.h"
|
|
#include "progressBar_imp.h"
|
|
// 命令规则 n为阶数 m为次数
|
|
class sph2xyz
|
|
{
|
|
public:
|
|
sph2xyz(){
|
|
Anm = NULL;
|
|
Bnm = NULL;
|
|
Pnm = NULL;
|
|
mCos = NULL;
|
|
mSin = NULL;
|
|
coff_C = NULL;
|
|
coff_S = NULL;
|
|
multi_array = NULL;
|
|
obs_topo = NULL;
|
|
}
|
|
~sph2xyz(){
|
|
if (Anm != NULL)
|
|
{
|
|
for (int i = 0; i < NN_size; i++)
|
|
{
|
|
if(Anm[i] != NULL) delete[] Anm[i];
|
|
if(Bnm[i] != NULL) delete[] Bnm[i];
|
|
if(Pnm[i] != NULL) delete[] Pnm[i];
|
|
if(coff_C[i] != NULL) delete[] coff_C[i];
|
|
if(coff_S[i] != NULL) delete[] coff_S[i];
|
|
}
|
|
delete[] Anm; delete[] Bnm; delete[] Pnm; delete[] coff_C; delete[] coff_S;
|
|
}
|
|
if (mCos != NULL)
|
|
{
|
|
for (int i = 0; i < lon_size; i++)
|
|
{
|
|
if(mCos[i] != NULL) delete[] mCos[i];
|
|
if(mSin[i] != NULL) delete[] mSin[i];
|
|
}
|
|
delete[] mCos; delete[] mSin;
|
|
}
|
|
if (multi_array != NULL)
|
|
{
|
|
for (int i = 0; i < lat_size; i++)
|
|
{
|
|
if(multi_array[i] != NULL) multi_array[i];
|
|
}
|
|
delete[] multi_array;
|
|
}
|
|
if (obs_topo != NULL) delete[] obs_topo;
|
|
}
|
|
int readSHC(char*,char*,char*); //读入球谐系数
|
|
int initMatrix(char*,char*,char*); //初始化相关的矩阵大小
|
|
int initObs(char*,char*,char*); //初始化观测点 如只有范围参数则只初始化经纬度位置
|
|
int calSolution(); //计算地形
|
|
void get_a_nm_array();
|
|
void get_b_nm_array();
|
|
void NALF_SFCM(double,double);
|
|
private:
|
|
double** Anm;
|
|
double** Bnm;
|
|
double** Pnm; //伴随勒让德函数系数 这个函数只和观测位置的纬度/余纬度相关 同一纬度只需要计算一次
|
|
double** mCos; //不同次数cos函数值 这个值只和观测位置的经度相关 行数为不同经度位置 列数为不同次数 矩阵维度即为经度个数*阶次 一般估算在1000*1000级别
|
|
double** mSin; //不同次数sin函数值 其他与上同
|
|
double** coff_S; //球谐系数sin参数
|
|
double** coff_C; //球谐系数cos参数
|
|
double** multi_array; //乘子矩阵
|
|
spoint* obs_topo; //计算地形是的观测位置 即计算半径值
|
|
|
|
double 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;
|
|
};
|
|
//读取球谐参数文件 文件名 起止阶次 列序列
|
|
int sph2xyz::readSHC(char* filename,char* para,char* orders)
|
|
{
|
|
ifstream infile;
|
|
infile.open(filename);
|
|
if(!infile)
|
|
{
|
|
cout << "can not open: " << filename << endl;
|
|
return 1;
|
|
}
|
|
|
|
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))
|
|
{
|
|
cout << "error ==> wrong parameter of " << para << endl;
|
|
return -1;
|
|
}
|
|
//识别列次序
|
|
int order[4];
|
|
if (4 != sscanf(orders,"%d,%d,%d,%d",&order[0],&order[1],&order[2],&order[3]))
|
|
{
|
|
cout << "error ==> wrong parameter of " << orders << endl;
|
|
return -1;
|
|
}
|
|
|
|
//按照最大阶数初始化下半三角矩阵
|
|
NN_size = n_end + 1;
|
|
coff_C = new double* [NN_size];
|
|
coff_S = new double* [NN_size];
|
|
for (int i = 0; i < NN_size; i++)
|
|
{
|
|
coff_C[i] = new double [i+1];
|
|
coff_S[i] = new double [i+1];
|
|
}
|
|
//初始化矩阵值为0
|
|
for (int i = 0; i < NN_size; i++)
|
|
{
|
|
for (int j = 0; j < i+1; j++)
|
|
{
|
|
coff_C[i][j] = coff_S[i][j] = 0.0;
|
|
}
|
|
}
|
|
|
|
int n,m; //行列号
|
|
double temp_d,temp_c,temp_s;
|
|
_1dArray temp_row;
|
|
string temp_str;
|
|
stringstream temp_ss;
|
|
while (getline(infile,temp_str))
|
|
{
|
|
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 = int(temp_row[order[0]]);
|
|
m = int(temp_row[order[1]]);
|
|
temp_c = temp_row[order[2]];
|
|
temp_s = 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 sph2xyz::initObs(char* r_para,char* i_para,char* refsys)
|
|
{
|
|
//解析经纬度范围
|
|
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))
|
|
{
|
|
cout << "error ==> wrong parameter of " << r_para << endl;
|
|
return -1;
|
|
}
|
|
else altitude = 0.0;
|
|
}
|
|
//解析间隔
|
|
if (2 != sscanf(i_para,"%lf/%lf",&dlon,&dlat))
|
|
{
|
|
cout << "error ==> wrong parameter of " << i_para << endl;
|
|
return -1;
|
|
}
|
|
//解析参考球
|
|
if (2 != sscanf(refsys,"%lf/%lf",&refr,&refR))
|
|
{
|
|
cout << "error ==> wrong parameter of " << refsys << endl;
|
|
return -1;
|
|
}
|
|
|
|
spoint temp_spoint;
|
|
double lon,lat;
|
|
lon_size = floor((lonmax-lonmin)/dlon) + 1;
|
|
lat_size = floor((latmax-latmin)/dlat) + 1;
|
|
obs_topo = new spoint [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.rad = refRadius(temp_spoint.lat,refr,refR) + altitude;
|
|
temp_spoint.val = 0.0;
|
|
obs_topo[i*lon_size+j] = temp_spoint;
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
//初始化矩阵
|
|
int sph2xyz::initMatrix(char* type,char* para,char* norType)
|
|
{
|
|
//初始化GM与R
|
|
if (strcmp(para,"NULL")) //如果para不为NULL则识别参数 否则将GM与R初始化为MAX_BDL
|
|
{
|
|
if (2 != sscanf(para,"%lf/%lf",&GM,&R))
|
|
{
|
|
cout << "error ==> wrong parameter of " << para << endl;
|
|
return -1;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
GM = R = MAX_BDL;
|
|
}
|
|
//初始化归一化类型
|
|
if (!strcmp(norType,"g")) norSum = 4.0*pi;
|
|
else if (!strcmp(norType,"m")) norSum = 1.0;
|
|
else
|
|
{
|
|
cout << "error ==> wrong parameter of " << norType << endl;
|
|
return -1;
|
|
}
|
|
//初始化伴随勒让德函数矩阵
|
|
Pnm = new double* [NN_size];
|
|
for (int i = 0; i < NN_size; i++)
|
|
Pnm[i] = new double [i+1];
|
|
//初始化sin和cos矩阵
|
|
mCos = new double* [lon_size];
|
|
mSin = new double* [lon_size];
|
|
for (int i = 0; i < lon_size; i++)
|
|
{
|
|
mCos[i] = new double [NN_size];
|
|
mSin[i] = new double [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*pi/180.0);
|
|
mSin[i][j] = sin(j*lon*pi/180.0);
|
|
}
|
|
}
|
|
//计算勒让德函数系数
|
|
get_a_nm_array();
|
|
get_b_nm_array();
|
|
//初始化乘子矩阵
|
|
multi_array = new double* [lat_size];
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
multi_array[i] = new double [NN_size];
|
|
}
|
|
//根据不同类型计算乘子参数和乘子矩阵
|
|
if(!strcmp(type,"t")) //topography
|
|
{
|
|
multi_factor = 1.0;
|
|
#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] = 1.0;
|
|
}
|
|
}
|
|
}
|
|
else if (!strcmp(type,"d")) //gravity disturbance
|
|
{
|
|
if (GM == MAX_BDL || R == MAX_BDL)
|
|
{
|
|
cout << "-g option must be set for gravitational calculation" << endl;
|
|
return -1;
|
|
}
|
|
|
|
multi_factor = 1e+5*GM/(R*R); // 将输出值转换为mGal
|
|
#pragma omp parallel for private(i,j) shared(R,lon_size) schedule(guided)
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
for (j = 0; j < NN_size; j++)
|
|
{
|
|
multi_array[i][j] = pow(R/obs_topo[i*lon_size].rad,j+2)*(j+1);
|
|
}
|
|
}
|
|
}
|
|
else if (!strcmp(type,"g")) //gravity anomaly
|
|
{
|
|
if (GM == MAX_BDL || R == MAX_BDL)
|
|
{
|
|
cout << "-g option must be set for gravitational calculation" << endl;
|
|
return -1;
|
|
}
|
|
|
|
multi_factor = 1e+5*GM/(R*R); // 将输出值转换为mGal
|
|
#pragma omp parallel for private(i,j) shared(R,lon_size) schedule(guided)
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
for (j = 0; j < NN_size; j++)
|
|
{
|
|
multi_array[i][j] = pow(R/obs_topo[i*lon_size].rad,j+2)*(j-1);
|
|
}
|
|
}
|
|
}
|
|
else if (!strcmp(type,"p")) //geo-potential
|
|
{
|
|
if (GM == MAX_BDL || R == MAX_BDL)
|
|
{
|
|
cout << "-g option must be set for gravitational calculation" << endl;
|
|
return -1;
|
|
}
|
|
|
|
multi_factor = 1e+5*GM/R; // 将输出值转换为mGal
|
|
#pragma omp parallel for private(i,j) shared(R,lon_size) schedule(guided)
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
for (j = 0; j < NN_size; j++)
|
|
{
|
|
multi_array[i][j] = pow(R/obs_topo[i*lon_size].rad,j+1);
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
cout << "error ==> unknown calculation type of " << type << endl;
|
|
return -1;
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
int sph2xyz::calSolution()
|
|
{
|
|
//计算
|
|
int i,j,n,m;
|
|
double temp_d,lat;
|
|
|
|
ProgressBar *bar = new ProgressBar(lat_size,"Process");
|
|
for (i = 0; i < lat_size; i++)
|
|
{
|
|
bar->Progressed(i);
|
|
lat = latmin + dlat*i;
|
|
//计算伴随勒让德函数 对于同一个纬度只需要计算一次
|
|
NALF_SFCM(90.0-lat,norSum);
|
|
//这里可以使用并行加速计算外层循环 内层计算因为是递归计算因此不能并行
|
|
//一种并行方案更快一些
|
|
#pragma omp parallel for private(j,n,m,temp_d) shared(i,multi_factor) 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]);
|
|
}
|
|
}
|
|
obs_topo[i*lon_size+j].val = multi_factor*temp_d;
|
|
}
|
|
}
|
|
|
|
cout << "# lon(deg) lat(deg) reference-radius(m) topography(m)|gravitational field(mGal)" << endl;
|
|
for (int i = 0; i < lon_size*lat_size; i++)
|
|
{
|
|
obs_topo[i].info();
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
void sph2xyz::get_a_nm_array()
|
|
{
|
|
int i,j;
|
|
Anm = new double* [NN_size];
|
|
for (i = 0; i < NN_size; i++)
|
|
Anm[i] = new double [i+1];
|
|
//向下列推计算
|
|
#pragma omp parallel for private(i,j) schedule(guided)
|
|
for (j = 0; j < NN_size; j++)
|
|
{
|
|
Anm[j][j] = 0; //对角线上的值直接给0 反正用不到
|
|
for (i = j+1; i < NN_size; i++)
|
|
{
|
|
Anm[i][j] = sqrt(((2.0*i-1)*(2.0*i+1))/((i-j)*(i+j)));
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
void sph2xyz::get_b_nm_array()
|
|
{
|
|
int i,j;
|
|
Bnm = new double* [NN_size];
|
|
for (i = 0; i < NN_size; i++)
|
|
Bnm[i] = new double [i+1];
|
|
//向下列推计算
|
|
#pragma omp parallel for private(i,j) schedule(guided)
|
|
for (j = 0; j < NN_size; j++)
|
|
{
|
|
Bnm[j][j] = 0; //对角线上的值直接给0 反正用不到
|
|
for (i = j+1; i < NN_size; i++)
|
|
{
|
|
Bnm[i][j] = sqrt(((2.0*i+1)*(i+j-1)*(i-j-1))/((i-j)*(i+j)*(2.0*i-3)));
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
void sph2xyz::NALF_SFCM(double theta,double norSum)
|
|
{
|
|
//赋初值给前两个对角线上的值
|
|
//norSum为1时第一个值为1/sqrt(4.0*pi),归一化值为1, norSum为4.0*pi时第一个值为4.0*pi/sqrt(4.0*pi)=1,归一化值为4.0*pi
|
|
Pnm[0][0] = sqrt(norSum)/sqrt(4.0*pi);
|
|
Pnm[1][1] = sqrt(3.0)*sin(theta*pi/180.0);
|
|
//计算对角线上的值 递归计算 不能并行
|
|
for (int i = 2; i < NN_size; i++)
|
|
Pnm[i][i] = sin(theta*pi/180.0)*sqrt(0.5*(2.0*i+1)/i)*Pnm[i-1][i-1];
|
|
|
|
//计算次对角线(m+1,m)上的值 递归计算 不能并行
|
|
for (int i = 0; i < NN_size-1; i++)
|
|
Pnm[i+1][i] = cos(theta*pi/180.0)*sqrt(2.0*i+3)*Pnm[i][i];
|
|
|
|
//声明系数和迭代变量
|
|
int i,j;
|
|
//这里可以使用并行加速计算外层循环 内层计算因为是递归计算因此不能并行
|
|
#pragma omp parallel for private(i,j) schedule(guided)
|
|
for (j = 0; j < NN_size-1; j++)
|
|
{
|
|
for (i = j+2; i < NN_size; i++)
|
|
{
|
|
Pnm[i][j] = Anm[i][j]*cos(theta*pi/180.0)*Pnm[i-1][j] - Bnm[i][j]*Pnm[i-2][j];
|
|
}
|
|
}
|
|
}
|
|
|
|
#endif |