gctl_tutorials/examples/gobser_tricone_ex.cpp

/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * 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 "gctl/core.h"
#include "gctl/io.h"
#include "gctl/potential.h"
#include "iostream"

int main(int argc, char *argv[])
{
	gctl::array<gctl::vertex3dc> msh_vert;
	gctl::array<gctl::grav_tri_cone> msh_cone;
	gctl::array<gctl::gravcone_para> cone_para;

	// read element from Gmsh file
	std::ifstream infile;
	gctl::open_infile(infile, "data/gobser_tricone/template", ".msh");
	gctl::read_gmsh_node(infile, msh_vert, gctl::Packed);
	gctl::read_gmsh_element(infile, msh_cone, msh_vert, gctl::Packed, nullptr);
	infile.close();

	// set origin for the cone
	gctl::vertex3dc ori_vert(gctl::point3dc(0.0, 0.0, 0.0), msh_vert.size());
	// create grav_tri_cone
	for (int i = 0; i < msh_cone.size(); i++)
	{
		msh_cone[i].set_origin(ori_vert);
	}

	//initialize the tensors
	gctl::callink_gravity_para(msh_cone, cone_para);

	// read topography data
	gctl::_2d_vector txt_content;
	gctl::array<gctl::point3ds> topo_sph;
	gctl::text_descriptor desc;
	gctl::read_text2vector2d("data/gobser_tricone/topo.txt", txt_content, desc);

	topo_sph.resize(txt_content.size());
	for (int i = 0; i < txt_content.size(); i++)
	{
		topo_sph[i].rad = 1e+5 + txt_content[i][2];
		topo_sph[i].lon = txt_content[i][0];
		topo_sph[i].lat = txt_content[i][1];
	}

	// initiate observation points
	double xmin = 30, xmax = 50, dx = 0.25;
	double ymin = 30, ymax = 50, dy = 0.25;
	int xnum = round((xmax - xmin)/dx) + 1;
	int ynum = round((ymax - ymin)/dy) + 1;
	gctl::array<gctl::point3ds> obs_ps(xnum*ynum);
	for (int j = 0; j < ynum; j++)
	{
		for (int i = 0; i < xnum; i++)
		{
			obs_ps[i + j*xnum].rad = 1e+5 + 200.0;
			obs_ps[i + j*xnum].lon = xmin + i*dx;
			obs_ps[i + j*xnum].lat = ymin + j*dy;
		}
	}

	// calculate base gravity
	gctl::array<double> out_obs;
	gctl::array<double> rho(msh_cone.size(), 1.0);
	gctl::gobser(out_obs, msh_cone, obs_ps, rho, gctl::Vz);

	// interpolate topography
	int m, n;
	double colat1, colat2, lon1, lon2;
	gctl::point3ds tmp_p;
	gctl::point3dc tmp_c;

	gctl::array<double> out_topo(msh_vert.size());
	for (int i = 0; i < msh_vert.size(); i++)
	{
		tmp_p = msh_vert[i].c2s();
		m = floor(tmp_p.lon - 28.0)/0.5;
		n = floor(tmp_p.lat - 28.0)/0.5;

		colat1 = 90.0 - (28.0 + n*0.5);
		colat2 = 90.0 - (28.0 + (n+1)*0.5);
		lon1 = 28.0 + m*0.5;
		lon2 = 28.0 + (m+1)*0.5;

		tmp_p.rad = gctl::sph_linear_interpolate_deg(colat1, colat2, lon1, lon2, 
			90.0-tmp_p.lat, tmp_p.lon, 
			topo_sph[m + 49*n].rad, 
			topo_sph[m + 1 + 49*n].rad, 
			topo_sph[m + 49*(n+1)].rad, 
			topo_sph[m + 1 + 49*(n+1)].rad);

		out_topo[i] = tmp_p.rad - 1e+5;

		tmp_c = tmp_p.s2c();
		msh_vert[i].set(tmp_c);
	}

	// recalculate tensors
	gctl::callink_gravity_para(msh_cone, cone_para);

	gctl::array<double> out_obs2;
	gctl::gobser(out_obs2, msh_cone, obs_ps, rho, gctl::Vz);

	for (int i = 0; i < out_obs.size(); i++)
	{
		out_obs2[i] -= out_obs[i];
	}

	for (int i = 0; i < out_obs.size(); i++)
	{
		std::cout << obs_ps[i].lon << " " << obs_ps[i].lat << " " << out_obs2[i] << std::endl;
	}
	return 0;
}