GEO_WU: 2013

2013年11月14日星期四

matlab 使用subplot时添加总title

使用subplot时，直接使用title命令时，只能将title添加在分图上，如果想给总图添加title应如何操作？

只需在最后一个subplot结束之后加上两行代码即可：
axes
set(gca,'visible','off','Title',text('string','\bfHello World!'))

2013年9月11日星期三

ubuntu12.04下安装matlab过程：

1. 下载matlab镜像文件（.iso），假设文件名为Matlab2013a.UNIX.iso
2. 挂载iso文件，命令为： sudo mount -o loop Matlab2013a.UNIX.iso /mnt
3. 进入挂载目录： cd /mnt
4. 安装: sudo ./install
5. 出现图形化安装界面，根据提示输入key以及license文件；
6. 默认安装位置为： /usr/local/MATLAB/R2013a
7. 更改权限： sudo chmod a+w -R ~/.matlab （重要）

8. 运行软件，命令为：sudo /usr/local/MATLAB/R2013a/bin/matlab （每次运行matalb程序输入命令太长，可以在.bashrc文件中加入alias matlab='bash /usr/local/MATLAB/R2013a/bin/matlab'，这样以后只要输入matlab命令即可打开软件）

注：由于license文件必须在ife的网络中才能有效，所以在以后的安装中注意这点，避免重复安装和卸载。

2013年5月30日星期四

自定义colormap meep (zt)

matlab 自定义colormap meep (from red through white to blue)

I have been thinking of using the colormap used in the MEEP software (an open source FDTD code) as a nice alternative to the classical Matlab colormap. Here is the short Matlab code to generate it and the resulting images:


mincolor    = [1 0 0]; % red
mediancolor = [1 1 1]; % white   
maxcolor    = [0 0 1]; % blue      

ColorMapSize = 4;
int1 = zeros(ColorMapSize,3); 
int2 = zeros(ColorMapSize,3);
for k=1:3
    int1(:,k) = linspace(mincolor(k), mediancolor(k), ColorMapSize);
    int2(:,k) = linspace(mediancolor(k), maxcolor(k), ColorMapSize);
end
meep = [int1(1:end-1,:); int2];

colormap(meep); 
colorbar

参考资料:

http://meyavuz.blogspot.de/2011/02/meep-colormap-for-matlab.html

2013年5月9日星期四

线性方程组求解与矩阵求逆

1. 线性方程组
在大地测量数据处理方面，线性方程组的求解无疑是最重要的问题之一。线性方程组的通用形式可写为$ AX=B $，系数矩阵A的维数为 $ (M，N) $，$ M $为观测方程的个数（等于观测值的个数，1个观测值一个方程），$N$为需要求解的系数个数，因此向量$ X $ 的长度为$ N $，观测值向量$ B $的长度为$ M $。根据$M$和$ N $的大小观测，可以将线性方程组分为三类：

$ M > N $, 超定方程组（也称矛盾方程组），此时通常求方程的LS解；
$ M = N $，适定方程组，有唯一解；
$ M < N $，欠定方程组，有无数组解，但通常求其最小范数解。

接下来主要讨论超定方程组和适定方程组的情况。

对于超定方程组，根据LS准则，可得到法方程：$ A^TAX=A^TB $，令$ C=A^TA，D=A^TB $，则$ C $的维数为$(N，N)$，D的维数为$(N，1)$，此时法方程的求解与适定方程组一致。下面的讨论以适定方程组为例。

2. 线性方程组的求解与求逆的关系

线性方程组$ AX=B $的求解在数学上表示为$ X=A^{-1} \times B $，但是在实际数值计算过程中，并不是一定要去计算$ A^{-1} $，如常用的消去法，三角分解前推回代法等。

矩阵求逆等价于求解方程组$ AX=I $，理论上与求解线性方程组一致。

在数值计算过程中，我们必须考虑计算的速度以及对内存的耗用情况，此时求逆和线性方程组的求解是两个不同的概念。

在求解线性方程组之前，必须明确你所需要的结果。如果只要最后结果$ X $，则尽量不要去求逆；

如果系数阵$ A $的求逆无法避免，尽量选择快速的求逆算法。当然，对于小矩阵，就不用考虑两者的差别了。

用matlab函数表示以上情况：

$ X=A\B $ (基于LU分解，前推回代直接求解参数向量$ X $)；

$ X=A^{-1} \cdot B $ （基于求逆求解参数向量$ X $ ）；

通常第一种方法速度快，占用内存少，但无法给出$ A^{-1} $的信息。

3. 线性方程组的常用数值解法

通常把线性方程组的数值解法分为直接法和间接法。所谓直接法，就是经有限次运算即可求得方程组的精确解的方法；迭代法将求解的方程组转化为一个无限序列，其极限就是方程组的解。直接法一般来说用于求解系数矩阵阶数不是太高的问题，工作量较小，精度较高，但程序较复杂；迭代法主要用于一些高阶问题，一般程序较为简单，但工作量有时较大，耗用内存较多。

1> 直接法

a. 克莱姆法则

$ x_i=\frac{det(A_i)}{detA} $

算法复杂度：$ (n+1) \times [(n-1) \times n! + n!] $

运算量非常大，理论上可行，实际计算中很少用；

三角形方程组的求解简单快速，所以通常将系数阵A转化为三角阵再进行运算。不同的分解形式对应不同的分解方法。求解下三角阵L的过程，称为前推$ (X_1-X_2---X_n) $;求解上三角阵U的过程，称为回代$ (X_n---X_2-X_1)$。前推和回代的计算量一样，为$ \frac{n^2}{2} $次乘法和加减法以及n次除法。

b. 高斯消去法

利用消去法，将方程组变换为上三角形方程组；整个消去过程需完成$ O(\frac{1}{3}n^3) $次乘法和加法运算；回代过程需$ \frac{n^2}{2} $次乘法和加减法以及n次除法，故用消去法求解n阶方程组的计算量约为$ O(\frac{2}{3}n^3+2n^2) $；

从高斯消去法进一步完善，产生了列主元和全主元消去法，这两种方法在消去过程中的计算量进一步增加，但方法上更严谨。

每一步消去过程对应系数阵A左乘了一个初等变换矩阵，这一系列初等变换矩阵的乘积对应一个下三角矩阵L，故高斯消去法理论上等价与LU分解法。

c. LU分解法

将系数阵A直接分解成两个三角阵的乘积，根据矩阵运算的规则计算每个三角阵的元素。分解后$ AX=B $的计算转化为两个三角形方程组的求解：

$ LY=B $ （前推，计算出Y）

$ UX=Y $ （回代，计算出X）

算法复杂度为 $ O(\frac{2}{3}n^3+2n^2) $，所以直接分解法与消去法的计算工作量基本上是相同的。

求解逆阵的算法复杂度为$ O(2n^3) $，是求解线性方程组的3倍。

d. 正交分解法

将矩阵A分解为一个正交阵和一个上三角阵的乘积，$ A=QR $。根据分解方法，通常有Householder方法（镜像映射法）和Schmidt正交化方法。

运算复杂度为：$ O(\frac{4}{3}n^3+2n^2) $，与上面两种方法相比，运算量多一倍。

该方法在求解超定方程组时，较为优越，原因如下：
$ A^TAX=A^TB $ with $ A=QR $, then $ R^TQ^TQRX=R^TQ^TB $ and finally, $ RX=Q^TB $，因此，QR分解大大减少了计算量。

e. 对称正定矩阵的平方根法和$ LDL^T$分解法
对称正定矩阵是特殊矩阵，可以用平方根法和$ LDL^T$分解法（通常也叫做Cholesky decomposition），其运算量和存储量较普通消去法和LU分解法均节省一半，为$ O(\frac{1}{3}n^3+2n^2) $。
平方根分解：$ A=LL^T $
$ LDL^T$分解：$ L $ 为单位下三角阵

2> 迭代法
对于阶数不高的线性方程组，直接法有效，但对于阶数很高的矩阵，有时需要迭代法来解决。迭代的规则不同，迭代法也就不同。通常分为：线性迭代和非线性迭代；一阶迭代和P阶迭代；定常迭代和非定常迭代。具体，可参考教科书，如“数值计算方法（冯康）”。
共轭梯度法：非线性迭代法，适用于对称正定阵，不需选取任何迭代参数，使用方便。
此外，对于高阶方程组的处理通常还采用分块矩阵处理方法。节省内存需要量，但增加了计算量，需要更多的计算时间（由于内外存交换数据需要时间）。

4. 矩阵求逆
除了必须要求矩阵的逆，否则尽量不采用求逆阵的方法来求解线性方程组。
求逆阵的算法复杂度为$ O(2n^3) $，大约是同等方法求线性方程组$ O(\frac{2}{3}n^3) $的3倍。
矩阵乘积的算法复杂度和矩阵求逆是等价的。
矩阵求逆当前最快方法的复杂度为$ O(n^{2.37}) $
高斯消去法，LU分解法，SVD等，复杂度都为$ O(n^3) $，差别在于前面的系数。

特殊矩阵的处理通常会有特殊的方法，能大大减少运算量和存储量。

2013年5月7日星期二

lapack 线性方程组最小二乘求解函数dgels

1. DGELS函数简介

DGELS函数是利用最小二乘求解线性方程组（超定或欠定）AX=B或A'X=B的函数，处于driver level。利用LU分解或者QR分解进行求解。

2. DGELS函数调用格式

DGELS(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)

输入参数：

TRANS（CHAR）：值为'N'，表示No transpose, 对应AX=B；值为'T',表示Transpose，对应A'X=B;

M,N（INT）: 系数矩阵A的维数，M>N,超定（下面只以超定方程组为例）,M<N，欠定；

NRHS（INT）：观测值向量B的列数（通常为1）；

A（Double）：系数阵；

LDA（INT）：同DGESV函数中的LDA，LDA>=max(1,M)，通常等于M；

B（Double）：观测值向量，当TRANS=‘N’时，B（M，NRHS），当TRANS=‘T’时，B（N，NRHS）；

LDB（INT）：同DGESV函数中的LDB，LDB>=max(1,M)，通常等于M；

LWORK（INT）：Length of Work, 数组WORK的长度，见http://www.netlib.org/lapack/lug/node118.html；

输出参数：

B（Double）：B(1:N)为计算的结果；B（N+1：M）误差信息，具体见函数说明；

A（Double）：M>N, QR分解矩阵信息；M<N, LU分解矩阵信息；

WORK（Double）：维数为LWORK，存储空间(WORKspace)信息，具体见说明http://www.netlib.org/lapack/lug/node120.html

INFO：=0，成功执行；=-i，第i个参数illegal；=i，A的第i个对角元素为0，A非满秩阵，无LS解。

3. DGELS函数举例：

http://www.nag.co.uk/lapack-ex/node45.html#DGELS

4. 同类型函数

DGELSY：利用完全正交分解计算AX=B的LS解
DGELSS：利用奇异值分解计算AX=B的LS解
DGELSD：利用奇异值分解（divide and conquer）计算AX=B的LS解

5. 参考资料

2013年5月6日星期一

lapack 线性方程组求解函数dgesv介绍

1. DGESV函数简介
DGESV函数是lapack函数库求解线性方程组AX=B的函数，处于drive level。利用LU分解进行求解，A=PLU，P为置换矩阵，L为单位下三角矩阵，U为上三角矩阵。

2. dgesv函数参数介绍
a) 形式
DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO)
输入参数：
N（INT）：线性方程组的数目，矩阵A的维数；
NRHS（INT）：矩阵B的列数（一般为1，单列）；
A（Double）：系数矩阵A，维数（LDA，N），通常LDA=N；
LDA（INT）：Leading Dimension of A, LDA>=max(1,N);
B（Double）：矩阵B，维数（LDB，NRHS），通常LDB=N，NRHS=1；
LDB（INT）： Leading Dimension of B，LDB>=max(1,N);
输出结果：
A（Double）：存储LU分解的L和U矩阵，L为单位下三角阵且对角线的单位值1没有存储；
IPIV（INT）：维数为N，利用该矩阵的值可以推求置换矩阵P，starting from a unit matrix E, E矩阵的第i行与E矩阵的第IPIV（i）行进行交换。
如IPIV=[2 2 3 4]，则P矩阵的结果为：

B（Double）：存储最终的求解结果X值；

INFO（INT）：求解状态参数。INFO=0，成功执行；INFO<0，则INFO=-i，第i个参数是illegal value；INFO>0, U(i,i)=0, LU分解成功执行，但是U是奇异的，最终结果无法计算；

DGESV函数调用了Lapack Computation level的两个函数DGETRF和DGETRS，前者是进行LU分解的函数，后者是利用分解的结果进行前向和后向推求最终结果X。简单介绍下这两个函数的调用格式：

DGETRF（M，N，A，LDA，IPIV，INFO），矩阵A的维度为（M，N），其它参数同上，返回结果为A（PLU），IPIV以及INFO。

DGETRS（TRANS，N，NRHS，A，LDA，IPIV，B，LDB，INFO），求解AX=B或A'X=B，该函数只有在INFO=0时才可以进行。输入变量'TRANS'为字符型，其意义如下：'N'对应AX=B； ‘T’对应A'X=B; 'C' 对应A'X=B (conjugate transpose)，其它变量的意义同DGESV函数。

3. dgesv函数举例
http://www.nag.co.uk/lapack-ex/node5.html
最终返回结果：
B=X= [ 1.0000 -1.0000 3.0000 -5.0000];
A=PLU=[
5.2500 -2.9500 -0.9500 -3.8000
0.3429 3.8914 2.3757 0.4129
0.3010 -0.4631 -1.5139 0.2948
-0.2114 -0.3299 0.0047 0.1314]

IPIV=[2 2 3 4]; 从一个单位阵开始，单位阵的第一行和第2行交换，第二行和第2行交换，第三行和第3行交换，第四行和第4行交换，

所以：

P=[0 1 0 0; 1 0 0 0; 0 0 1 0; 0 0 0 1];

L=[1 0 0 0; 0.3429 1 0 0; 0.3010 -0.4631 1 0; -0.2114 -.3299 0.0047 1];

U=[5.2500 -2.9500 -0.9500 -3.8000; 0 3.8914 2.3757 0.4129; 0 0 -1.5139 0.2948; 0 0 0 0.1314]

经过验证PLU=A；

4. dgesvx函数
dgesvx是expert dgesv函数，该函数可以返回更多的中间结果，如计算矩阵A的条件数，估计解的误差范围等等，相应的，函数也增加了很多参数，所以占用的内存也就更多。具体使用可以参见user book http://www.netlib.org/lapack/lug/node26.html和应用举例 http://www.nag.co.uk/lapack-ex/node6.html#DGESVX

2013年5月4日星期六

fortran编译器

Windows和Linux下Fortran编译器众多，这里简要介绍一些常用的编译器。

一. Windows平台下
在Win平台下，用的较多的是基于微软Visual Studio而集成的一整套IDE、商业编辑器：CVF和IVF。
1. CVF（Compaq Visual Fortran）
最常见的是CVF6.5，基于VS6，是win98时代的东西，更新到V6.6c被Intel收购了，即后来的IVF。CVF只支持到Fortran 95标准。
2. IVF（Intel Visual Fortran）
CVF的后继者，支持扩展的2003标准。IVF对Fortran77，90，95的兼容性不错。

二. Linux平台下
Linux下的fortran编译器很丰富，而且几乎都是免费的，如g77, gfortran, ifort 以及g95。
g77以前是GNU里的成员，后来被gfortran替代了。g95是在gfortan的发展过程中独立出去的，现在就只有作者一个人在维护。ifort是Intel的产品，对于单独的用户是免费的。这四个编译器用的较多的是gfortran和ifort，关于这两个编译器的讨论也集中在对老版fortran代码的兼容性（容错性）以及编译后的程序运行效率问题。总体来说，ifort的兼容性更好些，效率问题不太清晰。以下下面这两篇帖子对这个问题做过一些讨论，希望读完后自己能有一个总体的评判。

注：关于ifort在ubuntu下的安装可以参加另一篇博文：

http://geowu.blogspot.de/2012/07/ubuntuifort-xe.html

2013年4月18日星期四

latex宏包安装

在ubuntu下安装latex宏包过程如下：

下载宏包.sty文件或者是自己编写的.sty文件；
找到系统默认宏包的位置，通常是在/usr/share/texmf/tex/latex目录下；
将下载的宏包或自己编写的宏包移到系统对应的宏包位置；
更新一下： sudo texhash

举例：

自己编写的beamercolorthemeLUH.sty文件；
系统对应的beamercolortheme宏包位置为：/usr/share/texmf/tex/latex/beamer/base/themes/color
sudo mv beamercolorthemeLUH.sty /usr/share/texmf/tex/latex/beamer/base/themes/color
sudo texhash

2013年4月17日星期三

关于字符编码的一点知识

分享两篇介绍字符集和字符编码的博文：

此外，常用的latin字符集介绍见另一个链接：

https://zh.wikipedia.org/wiki/ISO/IEC_8859

2013年4月4日星期四

英文科技文献写作中的一些细节

在英文科技文献写作过程中，经常在一些细节上模棱两可，这里将这些细节做一汇总。

1. 参考文献和引用的撰写（基于journal paper）
不同的期刊有不同的参考文献格式，并没有一个统一的参考文献格式，这里只给出一种常用的格式。参考文献包含如下信息：作者+年份+标题+期刊名+卷号+页码，这里重点讨论作者的写法。首先列出常用的几条准则：

所有参考文献是按第一作者的last name (family name, or surname)的首字母顺序排列；
姓名缩写的规则是last name （用于正式场合）不缩写，given name （first name and middle name）可以缩写；
参考文献中必须包含所有在文章中出现的引用，and no more!
同一作者在一年里有多篇文章，需要在年份后加小写字母来区分这些文章；
直接引用别人的文章，必须加引用（但不推荐这种引用方法）；转述别人文章中的idea或者facts，也必须加引用（大多数为这种情况）。
不仅仅book和paper, internet source, computer software, e-mail and verbal communication都需要加以引用说明;
用了别人文章中的图，也必须加以引用；

接下来给出常用的参考文献和其对应的引用：

Single Author Bugjuice, B. 1970. Physiological ... (Bugjuice, 1970)

Two Authors Timm, T., Bugjuice, B. 1990. The role ... (Timm and Bugjuice, 1990)

Muitiple Authors Bugjuice, B., Timm, T., Cratchet, R. 1990. The ... (Bugjuice et at., 1990)

另外两种常见情况：

同时引用多个文章，如（Gumward, 1952; Bugjuice, 1970）。一个作者一年有多篇文章(Bugjuice, 1970a)

注：在写引用时，只给出last name就可以了。

参考资源：

http://abacus.bates.edu/~ganderso/biology/resources/writing/HTWcitations.html

http://tim.thorpeallen.net/Courses/Reference/Citations.html

2. 数字和单位之间是否要加空格？ Yes.

3. 单位符号是用斜体还是正体表示？正体。

4. 括号前后是否需要加空格？ Yes. 如： We spent NT$3,000 (about 100 U.S. dollars) on dinner. 但和括号里的第一个和最后一个字之间不加空格。

5. i.e., e.g. 和etc.的用法
i.e. 是拉丁文 id est 的缩写，它的意思就是“那就是说，换句话说”，等同于“that is，in other words” ，目的是用来进一步解释前面所说的观点。
e.g. 是拉丁文 exempli gratia 的缩写，它的意思是“举个例子，比如”，等同与“for example”,目的就是用几个例子来说明前面的观点。
etc.就比较好理解了，它是 etcetera 的缩写，意思是“等等”，相当于“and so on”
e.g. 和 etc. 不能出现在同一句话中，因为 e.g. 是表示泛泛的举几个例子，并没有囊括所有的实例，其中就已经包含“等等”，如果再加一个 etc. 就画蛇添足了，例如下面这句话就是错的：
Writing instructors focus on a number of complex skills that require extensive practice (e.g., organization, clear expression, logical thinking, etc.)
注意，用它们的时候，不要把“小点”给落了~

参考资源：http://bbs.ednchina.com/BLOG_ARTICLE_43177.HTM

2013年3月25日星期一

ubuntu下安装Lapack

1. 关于Lapack和Blas
Lapack (Linear Algebra Package)是一个以Fortran90编写，用于数值计算的函数集。例如，用于解线性方程组、计算特征向量、计算矩阵QR分解、奇异值分解等问题。
Blas(Basic Linear Algebra Subprograms)，基础线性代数程序集，用于基本的矢量和矩阵运算。最初发布于1979年，并用于建立更大的数值程序包（Lapack）。一级Blas负责矢量-矢量运算，二级Blas负责矢量-矩阵运算，而三级Blas负责矩阵-矩阵运算。由于Blas的高效，在高质量线性代数软件中常常使用它。

2. Lapack安装
1> 下载Lapack
http://www.netlib.org/lapack/ 最新版本3.4.2
2> 解压
tar -zxf lapack-3.4.2.tgz
3> 编辑make.inc文件
进入解压后的文件夹，会有一个make.inc.example文件，将该文件重命名为make.inc文件；并查看该文件，默认编译器为gfortran，根据需要修改对应的参数。在ubuntu下，采用gfortran作为编译器，一般不许要修改该文件。
4> 编辑Makefile文件
将下列语句（设定生成的库，默认不包含blaslib）
lib: lapacklib tmglib
# lib: blaslib variants lapacklib tmglib
改成：
# lib: lapacklib tmglib
lib: blaslib variants lapacklib tmglib
生成的库中包含了blaslib。
5> make
安装并测试Lapack。
6> 将生成的静态库(.a)移到/usr/lib文件夹下
生成的静态库包含libblas.a liblapack.a libtmglib等三个库。
sudo mv *.a /usr/lib

3. Lapack使用
在编译程序时，语句如下（其中test.f90文件中包含了lapack里的routines）：
gfortran -o test test.f90 -llapack -lblas

4. Lapack概述
1> Levels of Routines
Lapack的程序分三个层次，分别是driver， computational和auxiliary。driver routines解决一个具体的问题，如求解线性方程组或计算一个实对称阵的特征值。如果用户的问题刚好可以用该层次的程序解决，推荐用该层次的程序。这类程序通常调用几个computational routines或者一些其他的子程序。computational routines解决单独的一个计算任务，如LU分解。用户可以根据自身的需求调用该类程序。auxiliary routines这里不做阐述。一般用户一般调用driver和computational routines较多。
2> 数据类型和函数命名规则
Lapack函数命名格式为XYYZZZ，其中：
X表示数据类型：

S 表示单精度；
D 表示双精度；
C 表示复数；
Z 表示双精度复数；

YY表示矩阵类型：

DI 表示对角阵；
GE 表示一般矩阵；
GT 表示一般三角阵；
OR 表示正交阵；
PO 表示正定阵；
TR表示三角阵；

ZZZ表示具体的计算问题：

SV或SVX（求解线性方程组），前者是simple版，后者是expert版（输出一些中间结果，耗费内存更多）
QR（QR分解）
LU （LU分解）
...

3> 几个常用的函数名

dgesv 求解线性方程组AX=B
dgels 线性方程组AX=B的LS解
xyyTRF 将矩阵yy进行三角分解（根据矩阵yy的类型，选择不同的分解方法，如GE矩阵用LU分解，对称正定矩阵采用Cholesky分解）
xyyTRS 基于xyyTRF分解的结果求解AX=B
xyyTRI 基于xyyTRF分解的结果求解$A^{-1}$

因此，dgesv等价与dgetrf和dgetrs的组合，dposv等价与dpotrf和dpotrs的组合，dgels等价与dgeqrf和dtrtrs的组合。

更多内容参见Lapack User's Guide, http://www.netlib.org/lapack/lug/

5. 不错的资源

官方手册 http://www.netlib.org/lapack/lug/
用法举例 http://www.nag.co.uk/lapack-ex/lapack-ex.html#begin
中文帮助手册 http://wenku.baidu.com/view/d6c95d333968011ca30091fc.html
程序查询 http://www.netlib.org/lapack/explore-html/index.html

补充：
Lapack是用fortran写的，Lapacke是它的C语言接口，先安装Lapack,然后安装Lapacke,那么在程序中可以直接调用C函数来实现需要的功能。

2013年3月16日星期六

imagesc函数 NaN值颜色设定

imagesc函数主要是用来显示图片，其作用与pcolor类似，两者之间的比较见另一篇帖子。

这里主要介绍利用imagesc函数绘图时，矩阵值为NaN的通常用blue色来表示。如图1所示：

图1. NaN值对应为blue

如果想要将NaN值对应的颜色用白色表示，可以用如下语句：
h=imagesc(C);
set(h,'alphadata',~isnan(C));
所得新的图片如图2：

图2： NaN对应white

matlab图形边框不完整

在用matlab绘图时，有时会出现图形的边框缺少一边或两边，如图所示：

图1：图形边框不完整

解决此问题的方法如下：

方法一：将图片保存成jpg, png等格式时，边框可能会自动出现（不是确定管用的方法）；

方法二：设置图形的“渲染”属性，如：h=figure; set(h,'Renderer','opengl')；

方法三：开启opengl software，即添加语句：opengl software。

关于matlab函数opengl的用法，可以查看说明。

利用上述方法后，得到的图片如下：

图2 完整的图形边框

2013年2月23日星期六

格林尼治恒星时计算

参考以下几个网址：

注：通常还有一个IAU2000的GMST，其形式与常用的GMST略微有些不同，但计算结果相近。

2013年2月22日星期五

万有引力常量（GM）

常用的万有引力常量见如下网页：

2013年1月17日星期四

DE405星历计算天体位置和速度

按照如下步骤根据DE405天文星历计算天体的位置和速度：

1. 下载DE405星历数据和官方程序

下载地址：ftp://ssd.jpl.nasa.gov/pub/eph/planets/；

fortran文件夹里下载程序，ascii文件夹里下载星历数据。

2. 合并数据文件

在dos下： copy header.405+ascp2000.405+ascp2020.405 infile.405

在Unix下：cat header.405 ascp2000.405 ascp2020.405 > infile.405

3. 将合并后的数据文件转化为二进制文件

利用CV Fortran 编译asc2eph.f文件（由于源程序是用Fortran 77写的，以及平台差异，在Ubuntu编译或运行该程序时可以会有错误，但在CV Fortran下测试没有问题），生成asc2eph.exe文件,在dos下运行命令：asc2eph < infile.405. 运行的结果会生成二进制文件JPLEPH.

4. 测试源程序(testeph.f)

官方提供了testeph.f程序，测试源程序文件是否正确执行。在测试前需要将该程序文件中的部分参数赋值，具体如下：

在子程序FSIZER3中，将NRECL设为4，将NAMFIL设为JPLEPH，将KSIZE设为2036.在子程序STATE中，将语句CALL FSIZER3取消注释。保存修改后的程序。

继续在CV Fortran下编译程序testeph.f，生成可执行文件testeph.exe

进入dos下，执行命令：testeph < testpo.405

输出结果是一系列常用参数，以及和标准比较的结果，其中对应difference的位置数据全是E-13量级，此时表示程序运行正确；否则会出现"warning: next difference >= ..."字样，此时需要继续修改程序。

5. 将testeph.f文件中子程序“FSIZER3”以下所有的子程序拷贝到'selcon.f文件中'，并接下来准备在该文件前编写主函数。

6. 编写主函数

如果要计算星体的位置和速度主要是调用PLEPH函数：

PLEPH(JUD,NTARG,NCTR,R)，其中
输入参数为：JUD 儒略日； NTARG 目标星体的编号（参考函数说明）；NCTR 中心星体的编号；

返回参数为：R(6)，R(1:3)目标星体的位置（直角坐标），R(4:6)目标星体的速度。

默认结果的单位时AU.

注：由于编译环境的差别，源程序在Ubuntu下不一定能编译成功，但在CV Fortran下编译没有问题。

参考资料：
http://wenku.baidu.com/view/4cd0e24df7ec4afe04a1df13.html?from=related

GPS时转化为儒略日

GOCE PKI轨道的历元信息是用GPS时来表示的，如何将这些以秒为单位的GPS时转换成儒略日呢？以下将介绍一个简单的方法。

GPS时间系统的原点是1980年1月6日0时，对应的儒略日JUT0：

JUT0=2444244.5；

设GPS时为941068840.82s，则其相对GPS时间原点的间隔d为：

d=941068840.82/86400.0；

对应的JUT为：

JUT=JUT0+941068840.82/86400.

2013年1月11日星期五

fortran中字符和数字之间的相互转换

fortran中自带的ichar和char函数只能处理单个字符和数字之间的转化，在实际应用中经常会遇到多个字符和多位数字之间的转换，如：

字符串'500' 转换为整型500；
整型500转换为字符型'500'；
浮点型500.00转换为字符型；

1. 字符型转换为整型

例：'500'----->500

character(len=30) :: a

integer b

a='500'

read(a,*)b

write(*,'(I3)')b

2. 字符型转化为浮点型

例：'500.00'---->500.00

character(len=30) :: a

real b

a='500.00'

read(a,*)b

write(*,'(F6.2)')b

3. 整型转换为字符型

例：500----->'500'

character(len=3) :: a

integer b

b=500

write(a,'(I3)')b

write(*,*)a

4. 浮点型转换为字符型

例：500.13----->'500.13'

character(len=6) :: a