Posts

使用GMT程序包内的movie命令可以生成动画，更好地展示数据的动态变化。下面是一个利用该命令绘制动画的简单例子: 您的浏览器不支持视频标签 1. 准备程序 首先编写动画背景图的生成脚本，包括生成绘图所需的数据等工作也可以在此完成（如果已有外部数据则不需要）。在实际运行中不需要手动执行此脚本，而是通过movie命令调用。 # 使用heredoc方式将脚本保存到一个shell文件pre.sh # 在<<后添加-表示忽略tab制表符（注意不会忽略空格） cat <<- EOF > pre.sh # 使用gmt math生成数据，=号后接输出文件名称 # -T命令指定x坐标范围 # T SIND指定要计算的函数名称 gmt math -T0/360/10 T SIND = sin_point.txt gmt math -T0/360/1 T SIND = sin_curve.txt # 开始绘制底图 gmt begin # 使用gmt basemap绘制一张空的底图 具体的名称含义查看gmt文档 # -R指定坐标范围 -JX指定投影类型和地图大小 -X -Y平移图像 # -B指定坐标轴样式（包括ticks、labels和grids设置等） # --FONT_ANNOT_PRIMARY设置字体大小 gmt basemap -R0/360/-1.2/1.6 -JX22c/11.5c -X1c -Y1c \ -BWSne+glightskyblue -Bxa90g90f30+u@. -Bya0.5f0.1g1 --FONT_ANNOT_PRIMARY=9p gmt end # 结束 EOF 2. 准备主程序 编写动画生成的主脚本。movie命令内置了一批常量和变量可供使用，其中常量包括： ...

Introduction Data visualisation is extremely important for communicating the results of your research, either in a journal or to the general public, and for analysing and learning more about the characteristics of your data and system (so-called “exploratory data analysis”). One of the most fundamental tools in data visualisation is the two-dimensional plot (or graph). This tutorial will cover the basics of two-dimensional data visualisation using a program called gnuplot; a program which allows you to create high-quality, visually-pleasing figures and undertake robust post-hoc data analysis. ...

本系列共 4 篇：（一）Shell 技巧 · （二）怎样 getopts · （三）分享与扩展 · （四）C 语言 API 说明：本篇为归档笔记，内容仍在整理中，后续会继续补充。 简单的说明 作者在本系列的前几期教程中给大家分享了一些编写 GMT 脚本的小技巧，相信大家如果应用得当的话一定能有效提升你的脚本质量。本期教程所涉及的内容并不是 shell 脚本的相关内容，而是关于 GMT 动态库的使用方法。这也是我这几天研究的结果，想给大家粗略地介绍一下，也算给自己留个档。 本篇教程中所使用编程语言为 C/C++，因为也算是预编译脚本吧，姑且就放在这个系列里了。其实 GMT 的接口语言类型很丰富，包括 Matlab、Julia、python 等等，以后有机会再给大家分享这些语言的接口用法。目前，GMT 官方对于其动态库 API 的说明文档还不是很完善，因此很多地方我也是摸着石头过河，不是很清楚，请大家见谅。 本篇教程不涉及动态库的构建讲解，仅从一般使用方法入手。首先我们大致讲解一下 GMT 的 C 语言接口的总体使用策略。 总体策略 一般我们在命令行或者脚本中使用 GMT 命令时，典型的用法如下： $ gmt grd2cpt example.nc -Rd -Z -D > user.cpt 我们大致可以将这个命令分成这样几个部分： gmt grd2cpt 指定我们要使用的模块； example.nc 输入文件； -Rd -Z -D 模块使用的参数； user.cpt 输出文件。 其中，2 到 4 可以合并成一项，即某个模块运行需要的一系列参数，其实就是一个字符串。因此，我们首先祭出最核心的一个接口（函数）： int GMT_Call_Module (void *API, const char *module, int mode, void *args); 这个函数一共有 4 个参数，分别为： void *API GMT 模块调用总接口指针，同一个指针表示同一组调用。一般不用管，除非你在一个程序里要调用多组画图流程； const char *module 需要调用的模块名称，如 grd2cpt； int mode 调用的模式，一般为 GMT_MODULE_CMD（文后我们会列出可选的类型）； void *args 模块调用使用的参数，如 example.nc -Rd -Z -D > user.cpt。 看明白了么？GMT 的 C 语言接口其实非常的强大（傻瓜），我们要做的就是根据我们的数据编写相应模版的参数字符串。例如上面的 grd2cpt 命令从 API 调用就可以写成： ...

本系列共 4 篇：（一）Shell 技巧 · （二）怎样 getopts · （三）分享与扩展 · （四）C 语言 API 上一篇文章中我们介绍了如何使用 getopts 来获取 GMT 脚本执行时需要的命令行参数。通过使用 getopts 我们已经基本实现了 GMT 绘图过程的对象化，即给一个网格文件画一张图这样的基本操作。这样的工作流程的确很方便，但是它也有它不方便的地方。 其一是你不太可能在编写脚本之初就考虑到以后在使用过程中的所有需求，所以你可能还是需要经常为你的脚本添加和修改命令选项。这当然很好，你的脚本会越来越强大。但是过多的命令选项会让脚本的使用变得越来越困难，其他人可能需要花很多时间来学习它。而且，有的需求其实是一次性的，这时候去修改脚本显然并不划算。另一个不足之处是分享你的脚本将变得困难。设想这样一个情景，你想分享你的数据和对应的绘图脚本给你的同事，你发现你面对的并不是简单地将文件夹打包一下的情况。你需要将你的模版脚本拷贝到文件夹打包发送，然后告诉你的同事怎么使用这个脚本绘图，或者另写一个说明。这无疑增加了整个分享过程的复杂度。要是你的同事想将压缩包分享给第三人，那情况就更复杂了。 看到这里你可能会想“赶紧来点干的啊！”。 干的来了。 本篇教程我们就是要来解决分享与扩展性的问题，即如何将脚本特例化。什么意思？为了应对分享与扩展性这两个问题，我们的解决思路是在需要的时候为特定的数据（集）自动生成对应的绘图脚本。这个脚本内只包含 GMT 绘图命令，所以可以在任何装有 GMT 的电脑上运行。你的同事面对的只是单纯的 GMT 脚本，因此不需要额外学习任何其他知识。同时，你自己可以使用生成的脚本进行二次编辑，从而满足额外的绘图需求。这样就大大提高了模版脚本的扩展性。 说了这么多，到底怎么弄？ 答：写一个 Shell 函数，接收一个规则字符串和一个命令字符串，视规则执行或者打印命令字符串。就是说每当我们要执行一条 GMT 绘图命令时都要经过这个函数。如果我们想要绘图，这个时候就执行这个绘图命令；如果我们想要得到的是命令本身，就打印这条命令。 有了这个思路，我们来直接看这个函数： # 执行输入的语句或者将其显示在屏幕上 参数1（0显示 1执行）参数2（命令语句）参数3（命令语句） RunOrEcho() { # 预处理 去除命令语句中的制表符 first_str=`echo ${2// /''}` # 判断第三个参数是否存在 if [[ x${3} != x ]]; then sec_str=`echo ${3// /''}` fi # 如果第一个参数为1则执行命令语句 否则在屏幕上显示命令语句 if [[ ${1} == 1 ]]; then ${first_str} if [[ x${3} != x ]]; then ${sec_str} fi else if [[ x${3} != x ]]; then printf "%s\n%s\n" "${first_str}" "${sec_str}" else printf "%s\n" "${first_str}" fi fi } 简单地测试一下这个函数： ...

本系列共 4 篇：（一）Shell 技巧 · （二）怎样 getopts · （三）分享与扩展 · （四）C 语言 API 上一篇文章中我们介绍了一些在 GMT 脚本编写中可能会用到的 Shell 编程基础知识，今天我们将介绍 Linux 中一个很棒的内置命令 getopts。使用它我们可以方便地自定义标签并获取命令行参数。本篇教程将首先通过一个例子介绍 getopts 的使用方法，然后结合之前我们学习的 gmt_shell_functions.sh 脚本编写一个自动化绘制平面图的 Shell 脚本。如果你准备好了，我们就开始吧。 getopts 简短说明 话不多说，先上例子（personal_info.sh，不用着急看脚本，我们会在后面逐步讲解）： #!/bin/bash # 这个简单脚本会拾取姓名、年龄与性别三个参数并显示对应的信息。 # 初始化三个变量name age gender为Unknown name="Unknown" age="Unknown" gender="Unknown" # 显示个人信息 参数依次为 name age gender show_info() { if [[ ${2} -ge 18 && ${3} == 'm' ]]; then echo ${1} is a ${2} years ago man. elif [[ ${2} -ge 18 && ${3} == 'f' ]]; then echo ${1} is a ${2} years ago woman. elif [[ ${2} -lt 18 && ${3} == 'm' ]]; then echo ${1} is a ${2} years ago boy. elif [[ ${2} -lt 18 && ${3} == 'f' ]]; then echo ${1} is a ${2} years ago girl. else echo "Alien." fi } # 使用getopts从命令行获取参数 while getopts "hn:a:g:" arg do case $arg in # 拾取到-h选项，显示帮助信息后退出 h) printf "usage: ${0##*/} -n<name> -a<age> -g<gender> [-h]\n" printf "%s\t%s\n" "-n" "Your name" printf "%s\t%s\n" "-a" "Your age." printf "%s\t%s\n" "-g" "Your gender. Enter 'm' for male and 'f' for female." printf "%s\t%s\n" "-h" "Show help information." exit 0;; n) name=$OPTARG;; # 对变量进行赋值，注意后跟双分号结束。 a) age=$OPTARG;; g) gender=$OPTARG;; ?) # 拾取到未知参数，显示帮助信息后退出。 printf "error: unknow argument\nuse -h option to see help information.\n" exit 1;; esac done # 检查参数值并调用函数显示个人信息 if [[ ${name} != "Unknown" && ${age} != "Unknown" && ${gender} != "Unknown" ]]; then show_info ${name} ${age} ${gender} fi 在终端运行此脚本： ...

本系列共 4 篇：（一）Shell 技巧 · （二）怎样 getopts · （三）分享与扩展 · （四）C 语言 API 一般来说，我们在使用 GMT 绘制图件时可能会保留一些老的脚本。在需要时拷贝到工作文件夹中，再做适当修改并使用。这种方式的好处是对于特定数据有定制化的绘图脚本，可以方便的分享给他人。但是这种方式的不足即是脚本编写的工作量大、复用性差。而且随着脚本数量的增多，我们可能很难快速找到适用的老脚本。所以，利用 Shell 可编程属性，将常用的绘图流程标准化、模块化是提高 GMT 使用效率的重要一步，也是我们 GMT 脚本进阶的第一步。 本系列将分步骤介绍如何编写模块化的 GMT 脚本，帮助我们在日常使用中达到在任意文件夹内无需编写专门脚本也可以快速出图、批量化出图的要求。本篇内容为 Shell 编程中的一些基础知识，在我们后续的教程中会用到，大家可以先熟悉一下。 本教程中的脚本均在 Bash Shell 中测试通过。Windows 用户可使用 PowerShell。 使用变量 Shell 中可声明并使用变量，语法如下： # 注意变量名称与值之间用等号连接，不能有空格 name="Joe" age="19" # 使用变量时需添加$符号或使用${}包含，下面的语句将在终端显示Joe is 19 echo ${name} is ${age} # 对变量赋值时不需要添加$符号 age=20 # 下面的语句将在终端显示Joe is 20 echo ${name} is ${age} 为保证变量的正确调用，建议在任何时候都使用 ${ } 的方式。需要说明的是 Shell 中的所有变量都是作为字符串处理的，所以 age="19" 与 age=19 的含义是一样的。我们可以使用双引号或单引号来表示一个字符串，不同的是双引号中可包含变量，而单引号内的所有字符都将按照其原意进行处理。 使用函数 在 Shell 中我们可以将常用方法定义为函数，方便重复使用与管理。语法如下： # function标识可省略 function <func_name>() { <action> # return语句可省略 <return int> } 看一个具体的例子： ...

View on Github 线性共轭梯度算法库（C++ Library of Linear Conjugate Gradient，LIBLCG） 张壹（zhangyiss@icloud.com） 浙江大学地球科学学院·地球物理研究所 简介 liblcg是一个简单的C++线性共轭梯度算法库，包含了一般形式的共轭梯度与预优共轭梯度算法。可用于求解如下形式的线性方程组： Ax = B （可用共轭梯度求解）或者 PAx = PB（可用预优共轭梯度求解） 其中，A是一个N阶的对称矩阵、x为N*1的待求解的模型向量，B为N*1需拟合的目标向量。P则为预优算法中的预优矩阵，一般是一个N阶的对角阵。共轭梯度法广泛应用于无约束的线性最优化问题，拥有优良收敛与计算效率。 安装 GCC（MacOS） mkdir build cd build cmake .. make make install Windows 请自行拷贝代码新建项目并编译…嘿嘿😁！ 使用说明 自定义Ax计算函数 通常我们在使用共轭梯度法求解线性方程组Ax=B时A的维度可能会很大，直接储存A将消耗大量的内存空间，因此一般并不直接计算并储存A而是在需要的时候计算Ax的乘积。因此用户在使用liblcg时需要定义Ax的计算函数，同时提供初始解x与共轭梯度的B项（即拟合的对象）。如果使用预优方法还需要提供预优矩阵P项。Ax计算函数的形式必须满足算法库定义的一般形式： typedef void (*lcg_axfunc_ptr)(void* instance, const lcg_float* x, lcg_float* prod_Ax, const int n_size); 函数接受4个参数，分别为： void *instance 传入的实例对象； const lcg_float *input_array Ax计算中的x数组的指针； lcg_float *output_array Ax的乘积； const int n_size 矩阵的大小。 自定义进程监控函数 用户可以以下面的模版创建函数来显示共轭梯度迭代中的参数，并可以在适当的情况下停止迭代的进程。 typedef int (*lcg_progress_ptr)(void* instance, const lcg_float* m, const lcg_float converge, const lcg_para* param, const int n_size, const int k); 函数接收6个参数，分别为： ...

Program Propose 3D model construction and forward calculation of gravity and magnetic data using the Cartesian coordinates. density or magnetic models are built using elements of rectangular blocks. The program can either build 3D models or forward calculating gravity and magnetic data from a input model file. Some typical source types are supported by the program for fast model construction. File format of the 3D model used in this program is the 2.0 .msh file of the Gmsh software. ...

非线性共轭梯度法-DY法 算法可用于求解形如$min(f(x))$的非线性最优化问题，解的迭代形式为$x_{k+1}=x_k+a_kd_k$。其中$a_k$为迭代步长，$d_k$为迭代方向。（不同的$d_k$计算方法即为不同的非线性共轭梯度算法，它们对于$a_k$计算的要求也各有区别。） 要求目标函数$f(x)$连续可微且梯度$g(x)$已知。 DY法具有全局收敛性，但有一定几率陷入局部极小值。 DY法具有二次终止性，迭代终止于零梯度值位置或者$f(x) \leq threshold$处。 算法描述 step 1 -> 给出解的初值$f(x_0)$并计算其梯度$g(x_0)$，如果$\left\lVert g(x_0) \right\rVert = 0$则算法终止，说明$x_0$是$f(x)$的全局或局部极小解。否则设置$d_0 = -g_0$作为初始迭代方向（即解的最速下降方向）。 step 2 -> 执行满足如下所示的标准Wolfe条件的非精确线性搜索，计算得到$a_k$。 $f(x_k+a_k d_k)-f(x_k) \leq c_1 a_k {g_k}^T d_k$， ${g(x_k+a_k d_k)}^T d_K > c_2 {g_k}^T d_k$ 其中$0< c_1 < c_2 < 1$，一般取$c_1=0.1$，$c_2 \in (0.6,0.8)$。对于重磁非线形反演一般取$c_1=1e-4$，$c_2=0.9$。 step 3 -> 更新$x_{k+1}=x_k+a_kd_k $，若$\left\lVert g(x_{k+1}) \right\rVert = 0 $则算法终止，说明$x_{k+1}$是$f(x)$的全局或局部极小解。否则继续step 4。 step 4 -> 更新$d_{k+1}=-g_{k+1}+\beta_k d_k$，其中$\beta_k = \left\lVert g_{k+1} \right\rVert / {d_k}^T (g_{k+1}-g_k)$。$k:=k+1$返回step 2。 C++代码 func.h ...

View in toolbox 下载 stt 源码包 Program Propose The spherical triangular tessellation is a kind of partition of the spherical surface composed by only triangular cells. This program could generate the STT based on an icosahedron. The STT generated by the program could be refined around given points, lines, polygons and circles on the spherical surface. The exterior and interior outlines of the STT could also be customized. Installation You need to compile the executable file stt using the makefile. Please change the variable CC in the makefile to the compiler you use. ...