function [p,T,dT,Tm,e,ah,aw,la,undu] = gpt2_1w (dmjd,dlat,dlon,hell,nstat,it)
% (c) Department of Geodesy and Geoinformation, Vienna University of
% Technology, 2013
%
% The copyright in this document is vested in the Department of Geodesy and
% Geoinformation (GEO), Vienna University of Technology, Austria. This document
% may only be reproduced in whole or in part, or stored in a retrieval
% system, or transmitted in any form, or by any means electronic,
% mechanical, photocopying or otherwise, either with the prior permission
% of GEO or in accordance with the terms of ESTEC Contract No.
% 4000107329/12/NL/LvH.
% ---
%
% This subroutine determines pressure, temperature, temperature lapse rate,
% mean temperature of the water vapor, water vapor pressure, hydrostatic
% and wet mapping function coefficients ah and aw, water vapour decrease
% factor and geoid undulation for specific sites near the Earth surface.
% It is based on a 1 x 1 degree external grid file ('gpt2_1wA.grd') with mean
% values as well as sine and cosine amplitudes for the annual and
% semiannual variation of the coefficients.
%
% c Reference:
% J. B鰄m, G. M鰈ler, M. Schindelegger, G. Pain, R. Weber, Development of an
% improved blind model for slant delays in the troposphere (GPT2w),
% GPS Solutions, 2014, doi:10.1007/s10291-014-0403-7
%
% input parameters:
%
% dmjd: modified Julian date (scalar, only one epoch per call is possible)修改的儒略日日期(标量,每个调用只能有一个epoch)
% dlat: ellipsoidal latitude in radians [-pi/2:+pi/2] (vector)大地纬度的弧度
% dlon: longitude in radians [-pi:pi] or [0:2pi] (vector)经度的弧度
% hell: ellipsoidal height in m (vector)大地高
% nstat: number of stations in dlat, dlon, and hell同经纬度高度的站数
% maximum possible: not relevant for Matlab version
% it: case 1: no time variation but static quantities没有时间变化,只有静态量
% case 0: with time variation (annual and semiannual terms)随时间变化(年度和半年度时间段)
%
% output parameters:
%
% p: pressure in hPa (vector of length nstat) 气压hpa
% T: temperature in degrees Celsius (vector of length nstat)温度摄氏度
% dT: temperature lapse rate in degrees per km (vector of length nstat)温度递减率以每公里度为单位
% Tm: mean temperature of the water vapor in degrees Kelvin (vector of
% length nstat)水汽的平均温度
% e: water vapor pressure in hPa (vector of length nstat)水气压
% ah: hydrostatic mapping function coefficient at zero height (VMF1)
% (vector of length nstat)干映射函数
% aw: wet mapping function coefficient (VMF1) (vector of length
% nstat)湿映射函数
% la: water vapor decrease factor (vector of length nstat)水汽下降因子
% undu: geoid undulation in m (vector of length nstat)大地水准面起伏
%
% The hydrostatic mapping function coefficients have to be used with the
% height dependent Vienna Mapping Function 1 (vmf_ht.f) because the
% coefficients refer to zero height.
%
% Example 1 (Vienna, 2 August 2012, with time variation,grid file 'gpt2_1wA.grd):
%
% dmjd = 56141.d0
% dlat(1) = 48.20d0*pi/180.d0
% dlon(1) = 16.37d0*pi/180.d0
% hell(1) = 156.d0
% nstat = 1
% it = 0
%
% output:
% p = 1002.788 hPa
% T = 22.060 deg Celsius
% dT = -6.230 deg / km
% Tm = 281.304 K
% e = 16.742 hPa
% ah = 0.0012646
% aw = 0.0005752
% la = 2.6530
% undu = 45.76 m
%
% Example 2 (Vienna, 2 August 2012, without time variation, i.e. constant values):
%
% dmjd = 56141.d0
% dlat(1) = 48.20d0*pi/180.d0
% dlon(1) = 16.37d0*pi/180.d0
% hell(1) = 156.d0
% nstat = 1
% it = 1
%
% output:
% p = 1003.709 hPa
% T = 11.79 deg Celsius
% dT = -5.49 deg / km
% Tm = 273.22 K
% e = 10.26 hPa
% ah = 0.0012396
% aw = 0.0005753
% la = 2.6358
% undu = 45.76 m
%
%
% Klemens Lagler, 2 August 2012
% Johannes Boehm, 6 August 2012, revision
% Klemens Lagler, 21 August 2012, epoch change to January 1 2000
% Johannes Boehm, 23 August 2012, adding possibility to determine constant field
% Johannes Boehm, 27 December 2012, reference added
% Johannes Boehm, 10 January 2013, correction for dlat = -90 degrees
% (problem found by Changyong He)
% Johannes Boehm, 21 May 2013, bug with dmjd removed (input parameter dmjd was replaced
% unintentionally; problem found by Dennis Ferguson)
% Gregory Pain, 17 June 2013, adding water vapor decrease factor la
% Gregory Pain, 21 June 2013, using the 1 degree grid : better for calculating zenith wet delays (la)
% Gregory Pain, 01 July 2013, adding mean temperature of the water vapor Tm
% Gregory Pain, 30 July 2013, changing the method to calculate the water vapor partial pressure (e)
% Gregory Pain, 31 July 2013, correction for (dlat = -90 degrees, dlon = 360 degrees)
% Johannes Boehm, 27 December 2013, copyright notice added
% Johannes Boehm, 25 August 2014, default input file changed to
% gpt2_1wA.grd (slightly different humidity values)
% Johannes Boehm, 25 August 2014, reference changed to Boehm et al. in GPS
% Solutions
% ---
% change the reference epoch to January 1 2000
dmjd1 = dmjd-51544.5;
% mean gravity in m/s**2
gm = 9.80665;
% molar mass of dry air in kg/mol
dMtr = 28.965*10^-3;
% universal gas constant in J/K/mol
Rg = 8.3143;
% factors for amplitudes振幅的因子
if (it==1) % then constant parameters静态变量
cosfy = 0;
coshy = 0;
sinfy = 0;
sinhy = 0;
else
cosfy = cos(dmjd1/365.25*2*pi);
coshy = cos(dmjd1/365.25*4*pi);
sinfy = sin(dmjd1/365.25*2*pi);
sinhy = sin(dmjd1/365.25*4*pi);
end
% read gridfile读取格网数据
% use the 1 degree grid (GP)
fid = fopen('gpt2_1wA.grd','r');
% read first comment line读取第一行
line = fgetl(fid);
% initialization初始化
pgrid = zeros([64800, 5]);
Tgrid = zeros([64800, 5]);
Qgrid = zeros([64800, 5]);
dTgrid = zeros([64800, 5]);
u = zeros([64800, 1]);
Hs = zeros([64800, 1]);
ahgrid = zeros([64800, 5]);
awgrid = zeros([64800, 5]);
lagrid = zeros([64800, 5]);
Tmgrid = zeros([64800, 5]);
% initialization of new vectors初始化新向量
p = zeros([nstat, 1]);
T = zeros([nstat, 1]);
dT = zeros([nstat, 1]);
Tm = zeros([nstat, 1]);
e = zeros([nstat, 1]);
ah = zeros([nstat, 1]);
aw = zeros([nstat, 1]);
la = zeros([nstat, 1]);
undu = zeros([nstat, 1]);
% loop over grid points格网点的循环
% 64800 for the 1 degree grid (GP)
for n = 1:64800
% read data line
line = fgetl(fid);
vec = sscanf(line,'%f');
% read mean values and amplitudes
pgrid(n,1:5) = vec(3:7); % pressure in Pascal
Tgrid(n,1:5) = vec(8:12); % temperature in Kelvin
Qgrid(n,1:5) = vec(13:17)./1000; % specific humidity in kg/kg
dTgrid(n,1:5) = vec(18:22)./1000; % temperature lapse rate in Kelvin/m
u(n) = vec(23); % geoid undulation in m
Hs(n) = vec(24); % orthometric grid height in m
ahgrid(n,1:5) = vec(25:29)./1000; % hydrostatic mapping function coefficient, dimensionless
awgrid(n,1:5) = vec(30:34)./1000; % wet mapping function coefficient, dimensionless
lagrid(n,1:5) = vec(35:39); % water vapor decrease factor, dimensionless
Tmgrid(n,1:5) = vec(40:44); % mean temperature in Kelvin
end
fclose (fid);
% loop over stations站间的循环
for k = 1:nstat
% only positive longitude in degrees经度的度数只有正数
if dlon(k) < 0
plon = (dlon(k) + 2*pi)*180/pi;
else
plon = dlon(k)*180/pi;
end
% transform to polar distance in degrees转换成角度的极坐标距离
ppod = (-dlat(k) + pi/2)*180/pi;
% find the index (line in the grid file) of the nearest point查找最近点的索引(网格文件中的行)
% changed for the 1 degree grid (GP)
ipod = floor((ppod+1));
ilon = floor((plon+1));
% normalized (to one) differences, can be positive or negative规范化的差异,可�
评论7