initial upload

shc2xyz/CMakeLists.txt

@@ -0,0 +1,24 @@
+set(TOOL_NAME shc2xyz)
+set(BIN_DIR bin)
+set(INSTALL_DIR sbin)
+
+find_package(GCTL REQUIRED)
+include_directories(${GCTL_INC_DIR})
+
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O3")
+if(WIN32)
+	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3")
+else()
+	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3")
+endif()
+
+set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/${BIN_DIR})
+
+aux_source_directory(. TOOL_SRC)
+add_executable(${TOOL_NAME} ${TOOL_SRC})
+set_target_properties(${TOOL_NAME} PROPERTIES INSTALL_RPATH /usr/local/lib)
+set_target_properties(${TOOL_NAME} PROPERTIES CXX_STANDARD 17 CXX_STANDARD_REQUIRED ON)
+
+target_link_libraries(${TOOL_NAME} PUBLIC ${GCTL_LIB})
+
+install(TARGETS ${TOOL_NAME} RUNTIME DESTINATION ${INSTALL_DIR})

shc2xyz/func.h

@@ -0,0 +1,779 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * 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/algorithm.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

shc2xyz/main.cpp

@@ -0,0 +1,96 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * 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.
+ ******************************************************/
+
+#include "func.h"
+
+int main(int argc, char* argv[]) try
+{
+	gctl::flags_parser fp;
+	fp.set_proname("shc2xyz");
+	fp.set_proinfo("Forward modeling spherical harmonic coefficients to table data. \
+This program is a toolkit of the GCTL package. The GCTL comes with ABSOLUTE NO WARRANTY. \
+Please see instructions or contact the author for more information.");
+
+	fp.add_opt('t', "coeff-table", required_argument, NULL, "Input spherical harmonic coefficients.", "<filename>", true);
+	fp.add_opt('r', "range", required_argument, NULL, "Range of calculation.", "<lonmin>/<lonmax>/<latmin>/<latmax>[/<altitude>]", true);
+	fp.add_opt('i', "interval", required_argument, NULL, "Intervals of calculation.", "<dlon>/<dlat>", true);
+	fp.add_opt('d', "data-type", required_argument, NULL, "Objective data type. enter 'n' for null-specified (default), 't' for topography, 'd' for gravity disturbance,  'dlon' for gravity disturbance in the longitudinal direction, 'dlat' for gravity disturbance in the latitudinal direction, 'g' for gravity anomaly, 'p' for geo-potential, 'h' for height anomlay, and 'r' for radial gravity gradient", "<type>", true);
+	fp.add_opt('o', "order-degree", required_argument, NULL, "Start and end orders and degrees used for calculation.", "<ln>/<lm>/<hn>/<hm>", true);
+	fp.add_opt('s', "refer-system", required_argument, NULL, "Set reference system.", "<refr>/<refR>|WGS84|Earth|Moon|Mars", true);
+	fp.add_opt('l', "location", required_argument, NULL, "Read output locations from an input file, every line of the file has a longitudinal and a latitudinal value.", "<loc-file>", false);
+	fp.add_opt('g', "grav-para", required_argument, NULL, "GM and R multipliers for gravitational data's calculation.", "<GM>/<R>", false);
+	fp.add_opt('n', "normalization", required_argument, NULL, "Set normalized sum of associated Legendre functions. Enter 'm' for mathematic normalization which equal 1 and 'g' for geodetic normalization which equal 4*pi (default).", "<type>", false);
+	fp.add_opt('c', "columns", required_argument, NULL, "Select data columns of the input table data. the default is 0,1,2,3.", "<col1>,<col2>,<col3>,<col4>", false);
+	fp.add_opt('j', "jump-lines", required_argument, NULL, "Jump the first <num> lines of the input file.", "<num>", false);
+	fp.add_opt('a', "altitude", required_argument, NULL, "Provide a input file for observations' altitudes. use +d to select input data columns, the default is 0,1,2", "<alti-file>[+d<col1>,<col2>,<col3>]", false);
+	fp.add_opt('h', "help", no_argument, NULL, "Show help information.", 0, false);
+	fp.configure(argc, argv);
+
+	if(argc == 1 || fp.set_opt('h'))
+	{
+		fp.show_help_page(std::clog, " > <outfile>");
+		return 0;
+	}
+
+	std::string infilename, interfilename, range, interval;
+	std::string type, referSystem, coffPara, columns, jumps, gravPara, norType, altiFile;
+
+	fp.get_argv({'t', 'r', 'i', 'd', 'o', 's', 'l', 'g', 'n', 'c', 'j', 'a'}, 
+		{&infilename, &range, &interval, &type, &coffPara, &referSystem, &interfilename, &gravPara, &norType, &columns, &jumps, &altiFile});
+	if (type == "NULL") type = "n";
+	if (columns == "NULL") columns = "0,1,2,3";
+	if (norType == "NULL") norType = "g";
+	if (jumps == "NULL") jumps = "0";
+
+	// 查看是否通过强制参数检查
+	if (!fp.pass_mandatory()) return 0;
+
+	shc2xyz sx;
+	//读入球谐系数文件
+	if(sx.readSHC(infilename.c_str(), coffPara.c_str(), columns.c_str(), jumps)) return 0;
+	//初始化观测点位置
+	if(sx.initObs(range.c_str(), interval.c_str(), referSystem.c_str())) return 0;
+	//重定位观测点高程
+	if(sx.relocateAltitude(altiFile.c_str())) return 0;
+	//确定计算类型
+	if(sx.initTargetType(type.c_str())) return 0;
+	//初始化矩阵
+	sx.initMatrix(type.c_str(), gravPara.c_str(), norType.c_str(), altiFile.c_str());
+	//计算
+	if (altiFile == "NULL")
+	{
+		sx.calSolution();
+	}
+	else sx.calSolution2(type.c_str());
+	//输出
+	sx.outObs(interfilename.c_str());
+	return 0;
+}
+catch(std::exception &e)
+{
+	GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}