รองศาสตราจารย์ ดร.บุญทรัพย์ วิชญางกูร คณะวิศวกรรมศาสตร์ มหาวิทยาลัยธรรมศาสตร์

รองศาสตราจารย์ ดร.บุญทรัพย์ วิชญางกูร
Boonsap Witchayangkoon
(PhD, Associate Professor)

Engineering Building, Thammasat University
Pathumtani, THAILAND
Tel : (66)02-564-3001-9 ext. 3101, 3039
wboon @ engr.tu.ac.th

Education

• PhD (Spatial Information Science and Engineering) 2000, University of Maine, USA
• B.Eng. (Honors) (Civil Engineering) 1992, King Mongkut's Institute (now University) of Technology Thonburi, Thailand
• LL.B. 2007, Sukhothai Thammathirat Open University.

Committee Member (partial list):

• Committee Member of the Thammasat University Board of Trusstee 2005-2007 กรรมการสภามหาวิทยาลัยธรรมศาสตร์ 2548-50
• Committee Member of the Faculty Senate of Thammasat Univ 2002-now. กรรมการสภาอาจารย์มหาวิทยาลัยธรรมศาสตร์ 2545-ปัจจุบัน

Scholar Visitor

• Visiting Faculty (UMAP) at the Geoscience Institute, Shizuoka University, Japan: March 15- June 10, 2006

Fields of Interest
Spatial Information Engineering - วิศวกรรมภูมิสารสนเทศ
Informational Construction Management - สารสนเทศในการจัดการงานก่อสร้าง
GPS/GIS & analysis -
Transportation & Information of safety & accidents analysis
Positioning of military weapons & cannon explosions
Operations research in civil engineering
Basic Research in Civil Engineering & related fields
GPS-Aided Air Traffic Management
GPS application PPP for studying of Troposphere and Ionopshere
Engineering & Law - law-aided technology, environmental law, ..

On-Going Research Work

• GPS Supplementary Artificial Rainmaking Via Precipitable Water Vapor Information from Tropospheric Analysis


• Traffic Calming System on Thammasat University - e.g. Analysis of Speed Bump & ..

Graduate Students Research Group

• สัญญา นามี Sanya Namee, PhD Candidate


• อภิรัฐ ฉัตรานุสรณ์ Master Candidate

Publications

• Namee, S. and Witchayangkoon, B. 2011. Crossroads Vertical Speed Control Devices: Suggestion from Observation. International Transaction J of Engineering, Management, and Applied Sciences and Technologies, 2(2), 161-171.
  
• Sirimontree, S., B.Witchayangkoon, N. Khaosri, and J.Teerawong. 2011. Shear Strength of Reinforced Concrete Bream Strengthened by Transverse External Post-tension. American J of Engineering and Applied Science, 4(1), 108-115.
  
• Witchayangkoon, B., S.Namee, V.Temrungsri, and S.Intaratad. 2011. Traffic Calming at Crossroads Using Speed Bumps/Humps: A Case Study in Thammasat University, Rangsit Campus, Thailand. J of Science and Technology, 19(1) (Jan-Mar), 60-71. (in Thai language)
  
• Witchayangkoon, B. and P. Dumrongchai. 2009. Computational Linearized Least Square for 2D Positioning Analysis Using MathCad Programming. KMUTT Research and Development Journal. (Accepted)
  
• Sritipun, P. and B. Witchayangkoon. 2009. GIS Database of Wind Speed for the central region of Thailand. Proc. the 14th National Convention on Civil Engineering, Section of Survey and Geographic Information System Engineering (SGI), Nakorn-Ratchasrima, THAILAND, May 13-15, 2009. (Accepted)
  
• Namee, S., and B. Witchayangkoon. 2008. An Application of Geographic Information System for Assessment Risk Area from Hazardous Materials Transport. Proc. the 5th National Transport Conference (section Transport System & Sustainable Development Strategies), Thailand, Dec 19, 2008: 76-83.
  
• Witchayangkoon, B. 2004. A GIS Database of Pile Dynamic Load Tests for the Bangkok Region. Proc. International Conference on Civil and Environmental Engineering (ICCEE), Higashi-Hiroshima, Japan.
  
• Witchayangkoon, B. 2003. Estimation of Zenith Troposphere Using Dual-Frequency GPS Precise Point Position. Proc. the 4th Regional Symposium on Infrastructure Development in Civil Enigineering, Bangkok, Thailand.
  
• Wattanakul, N. and B.Witchayangkoon. 2002. "A survey and analysis of energy consumption of industrial factories within the central region of Thailand." Research and Development Journal of the Engineering Institute of Thailand, Vol. 13(4): 11-15.
  
• Witchayangkoon, B. "GPS-based vehicular velocity determination on the Chalermmahanakhon Expressway." Proc. the First Conference on Civil and Environmental Engineering (ICCEE-2002). Higashi-Hiroshima, Japan. Oct 30-31, 2002: 215-221.
  
• Witchayangkoon, B. 2002. Roles of GIS-aided construction management. Proc. the 8th National Convention on Civil Engineering. Hotel Sofitel Raja Orchid Khon Kaen, Thailand. Oct 23-25, 2002: CEM129-CEM134.
  
• Tounnawarat, N., B.Witchayangkoon, and K. Chompooming. 2002. "GPS/GIS application for historical monuments around the Rattanagosin Island." Proc. the 8th National Convention on Civil Engineering. Hotel Sofitel Raja Orchid Khon Kaen, Thailand. Oct 23-25.
  

Courses Taught @ Faculty of Engineering @ Thammasat University
Surveying
Engineering Mechanics
Strength of Material
Geographic Information System (GIS)
Global Positioning System (GPS)
Research Methodology
Applied information system
Spatial information analysis

Professional Services

• Teachings and research advising for national & international institutes and universities
• Consultancy assistance for many academics and professional projects.
• Devoted advisor for students projects for building or repairing schools in the very rural areas of Thailand since 2001.
  
• Professional Engineer (PE)
• Contributing reader of various national and international journals
• Energy & Environmental Conservation
• อาจารย์หลักสูตร วิศวกรรมศาสตร์สองสถาบัน (TEP) คณะวิศวกรรมศาสตร์ ม.ธรรมศาสตร์

This is a mirror homepage of Dr.Boonsap Witchayangkoon official web site @ Faculty of Engineering Thammasat University @ THAILAND

เวปนี้ เป็นเวปสำรองของ รองศาสตราจารย์ ดร.บุญทรัพย์ วิชญางกูร ภาควิชาวิศวกรรมโยธา คณะวิศวกรรมศาสตร์ มหาวิทยาลัยธรรมศาสตร์
http://www.ce.engr.tu.ac.th/ce/eng/staff_view.php?company_code=145


THAILAND : We are here..
For the popular international engineering education, we offer the attractive
Twining Engineering Program (TEP) and Thammasat English Program of Engineering (TEPE). Our programs are worldwide recognized to international standard of top-rank Thailand engineering education.

ELEMENTS OF GPS PRECISE POINT POSITIONING

ELEMENTS OF GPS PRECISE POINT POSITIONING
PhD Dissertation By Boonsap Witchayangkoon.
PhD Thesis Adviser: Professor Dr. Alfred Leick.
You can download freely from
http://www.cctechnol.com/uploads/167_0.pdf

*** My dissertation work has been mentioned in the book
GPS Satellite Surveying (see page 254-256).

Extended Kalman Filter for Analyzing GPS GNSS PPP

function [Result,Residual,Relativity,Pminus,CLKphi,NumSV,W,Pkalman] ...
= EKF(DOY,rTime,Ph1,Ph2,Pr1,Pr2,SP3X,SP3Y,SP3Z,SP3clk,Vx,Vy,Vz,x,Elev,trop,type)

% Run Extend Kalman Filtering
%
% %%%%%%%%% HELP %%%%%%%%%%
%
% [Result,Residual,Relativity,Pminus,CLKphi,NumSV,W,Pkalman] ...
% = EKF(DOY,rTime,Ph1,Ph2,Pr1,Pr2,SP3X,SP3Y,SP3Z,SP3clk,Vx,Vy,Vz,x,Elev,trop,type)
%
% INPUT
%
% DOY = Day of year [1 x 1]
% rTime = GPS second of RINEX file [GPSsec], [numep x 1]
% Ph1 = Carrier phase L1,cycles [numep x numsat]
% Ph2 = Carrier phase L2,cycles [numep x numsat]
% Pr1 = Pseudorange C/A code, metres [numep x numsat]
% Pr2 = Pseudorange P code , metres [numep x numsat]
% SP3X, SP3Y, SP3Z = Satellite coordinate (WGS84), metres [numep x numsat]
% SP3clk = Precise clock correction, microseconds [numep x numsat]
% Vx = Satellite velocity [numep x numsat], m/s
% Vy = Satellite velocity [numep x numsat], m/s
% Vz = Satellite velocity [numep x numsat], m/s
% x = Aproximate station coordinates, metres [3 x 1]
% Elev = Elevation angle, radians [numep x numsat]
% trop = 1 Estimate tropospheric correction from EKF
% = 2 Approximate tropospheric correction from Saastamoinen model
% = 3 Approximate tropospheric correction from Saastamoinen model & Neill mapping function
% type = 1 Normal
% = 2 transition matrix (Phi) is identity
%
% OUTPUT
%
% Result = [x y z Ni dT Tropo PDOP TDOP GDOP]
% Residual
% Relativity
% Pminus
% CLKphi
% NumSV
% W
% Pkalman
%
% Modified by Phakphong Homniam
% December 25, 2002
% Orginal source code in Mathcad by Boonsap Witchayangkoon, 2000

%%%%%%%%%% CONSTANT DECLARATION %%%%%%%%%%

% Global constants of the relevant parameters in constant_global.m file
constant_global

global NumSat NumPara gamma maskangle

% Satellite mask angle
maskangle = 15*deg;

% Kalman initialization constants
SDx = 0.0003; % (m) used in Q (dynamic model)
SDclk = 30; % (m) usde in Q (dynamic model)
SDph = 0.02; % (m) STDS of carrier phase (R)
SDpr = 0.2; % (m) STD of pseudorange (R)
SDcoorinit = 1; % (m) STD of initial parameters (P)
SDambinit = 10^10; % (m) STD of initial ambiguities (P)
SDclkinit = 500; % (m) STD for initial clock (P)
SDtropinit = 0.1; % (m) STD for initial tropo (P)

% Specification when estimating troposphere
% Tropospheric Correlation Time
thor = 18000; % seconds

% Number of satellite
NumSat = size(Ph1,2);

gamma = (f1/f2)^2;
tep = length(rTime);
TimeInterval = rTime(3)-rTime(2);

%%%%%%%%%% BEGIN %%%%%%%%%%

Result = [];
Residual = [];
Relativity = [];
Pminus = [];
CLKphi = [];
NumSV = [];
W = [];
Pkalman = [];

% Estimate receiver clock error for the first epoch for EKF
for i = 1:NumSat
if (Ph1(1,i) ~= NA) & (Pr1(1,i) ~= NA)
if SP3clk(1,i) ~= NAJ
rho = sqrt((SP3X(1,i)-x(1))^2+(SP3Y(1,i)-x(2))^2 ...
+(SP3Z(1,i)-x(3))^2);
h = rho-c*SP3clk(1,i)*10^-6;
FirstClkEst = Pr1(1,i)-h;
break
end
end
end

switch trop
case 1
% Number of unknown parameters
% (x,y,z)+(ambiguities)+Receiver Clock+Tropo
NumPara = 3+NumSat+2;

% Set martix Xm
zenith = 2.3;
Xstart = [x(:);zeros(NumSat,1);FirstClkEst;zenith];
Xm = Xstart;
Xminus = Xm;

% Set matrix Q
qSD(1:3) = SDx;
qSD(NumSat+4) = SDclk;
qSD(NumSat+5) = 0.01/3600; % may be 0.1/60 or 0.01/3600 : i don't sure
Q = diag(qSD.^2);

% Set matrix Pminus
pmSD(1:3) = SDcoorinit;
pmSD(4:NumSat+3) = SDambinit;
pmSD(NumSat+4) = SDclkinit;
pmSD(NumSat+5) = exp(-TimeInterval/thor);
Pm = diag(pmSD.^2);
Pminus = Pm;

case {2,3}
% Number of unknown parameters
% (x,y,z)+(ambiguities)+Receiver Clock
NumPara = 3+NumSat+1;

% Set martix Xm
Xstart = [x(:);zeros(NumSat,1);FirstClkEst];
Xm = Xstart;
Xminus = Xm;

% Set matrix Q
qSD(1:3) = SDx;
qSD(NumSat+4) = SDclk;
Q = diag(qSD.^2);

% Set matrix Pminus
pmSD(1:3) = SDcoorinit;
pmSD(4:NumSat+3) = SDambinit;
pmSD(NumSat+4) = SDclkinit;
Pm = diag(pmSD.^2);
Pminus = Pm;
end

% Set matrix Phi
switch type
case 1
Phi = eye(NumPara);
% Set Phi for the receiver clock
% = 0 pure white noise
% = 1 random walk
Phi(NumSat+4,NumSat+4) = 0;
case 2
Phi = eye(NumPara);
end
[PR,Phase] = typeobs(Ph1,Ph2,Pr1,Pr2,tep);
ep = 1;

for i = 1:(tep-1)
phase = Phase(i,:);
pr = PR(i,:); prp1 = PR(i+1,:);
sp3x = SP3X(i,:); sp3xp1 = SP3X(i+1,:); vx = Vx(i,:);
sp3y = SP3Y(i,:); sp3yp1 = SP3Y(i+1,:); vy = Vy(i,:);
sp3z = SP3Z(i,:); sp3zp1 = SP3Z(i+1,:); vz = Vz(i,:);
sp3clk = SP3clk(i,:); sp3clkp1 = SP3clk(i+1,:); elev = Elev(i,:);

% set matrix H
[HhZ1,HhZ2,HhZ3,HhZ4] = MatrixHhZ(Xm,phase,pr,sp3x,sp3y,sp3z,sp3clk,elev,x,DOY,vx,vy,vz,trop);
H = HhZ1;

% Set matrix R
rSD(1:NumSat) = SDph;
rSD(NumSat+1:2*NumSat) = SDpr;
R = diag(rSD.^2); % Observation eorror
% Reweight observation at low elevation angles
% for j = 1:NumSat
% if (elev(j) ~= NA) & (elev(j) < 30*deg)
% if elev(j) < 25*deg
% lv = elev(j)+5*deg;
% else
% lv = elev(j);
% end
% R(j,j) = (SDph/sin(lv))^2;
% R(j+NumSat,j+NumSat) = (SDpr/sin(lv))^2;
% end
% end

K = Pm*H'*inv(H*Pm*H'+R); % Calculate Kalman Gain
w = K*HhZ2;
W = [W,w];
Xestimate = Xm+w; % Using measurement to update estimate
P = Pm-K*H*Pm; % Update the error covariance
Pkalman = [Pkalman;P];
% Phi clk here
Phi(NumSat+4,NumSat+4) ... % Xestimate or x???
= clkphiall(Xestimate,pr,sp3x,sp3y,sp3z,sp3clk,prp1,sp3xp1,sp3yp1,sp3zp1,sp3clkp1);
Xm = Phi*Xestimate; % Project the state ahead
Xminus = [Xminus,Xm];

switch type
case 1
Pm = PhiPPhiT(Phi,P,trop)+Q; % Project the covariance ahead
case 2
Pm = P+Q;
end

Pminus = [Pminus;Pm];
Result(1:NumPara,ep) = Xestimate;
A = H4DOPs(phase,pr,sp3x,sp3y,sp3z,elev,x);
ATA = A'*A;
detATA = det(ATA);
if detATA == 0
PDOP = NA;
TDOP = NA;
GDOP = NA;
else
Pvar = inv(ATA);
PDOP = sqrt(Pvar(1,1)+Pvar(2,2)+Pvar(3,3));
TDOP = sqrt(Pvar(4,4));
GDOP = sqrt(Pvar(1,1)+Pvar(2,2)+Pvar(3,3)+Pvar(4,4));
end
Result(NumPara+1,ep) = PDOP;
Result(NumPara+2,ep) = TDOP;
Result(NumPara+3,ep) = GDOP;
Residual = [Residual,HhZ2];
NumSV(ep) = HhZ3;
Relativity = [Relativity,HhZ4(:)];
CLKphi(ep) = Phi(NumSat+4,NumSat+4);
ep = ep+1;
end

%%%%%%%%%% RESULT %%%%%%%%%%

Result = Result';
Residual = Residual';
Relativity = Relativity';
Pminus;
CLKphi = CLKphi(:);
NumSV = NumSV(:);
W = W';
Pkalman;
%%%%%%%%%% END %%%%%%%%%%

Neill Mapping Function Source Code for GPS/ GPS PPP/ GNSS PPP

function hmf = NeillMF(DOY,lat,h,elev)

% Neill Mapping Function
%
% %%%%%%%%%% HELP %%%%%%%%%%
%
%hmf = NeillMF(DOY,lat,h,elev)
%
% Input Data
% DOY is day of year.
% lat is latitude, degrees.
% long is longitude, degrees.
% h is height, meters.
% elev is elevation angle, degrees.
%
% Written by Phakphong Homniam
% September 20, 2002
% Original Mathcad source code by Boonsap Witchayangkoon, 2000

% Coefficient for changing DOY to radian
day2rad = 2*pi/365.25;

% Coefficient for changing degrees to radians
deg2rad = pi/180;

% Latitude array for the hydrostatic mapping function coefficients
lat_hmf = [15,30,45,60,75]';

% Average of coefficients a, b and c corresponding to the given latitude.
% See Table 5.4
abc_avg = [1.2769934 2.9153695 62.610505
1.2683230 2.9152299 62.837393
1.2465397 2.9288445 63.721774
1.2196049 2.9022565 63.824265
1.2045996 2.9024912 64.258455]*10^-3;

% Amplitude of coefficients a, b and c corresponding to the given latitude.
abc_amp = [0 0 0
1.2709626 2.1414979 9.0128400
2.6523662 3.0160779 4.3497037
3.4000452 7.2562722 84.795348
4.1202191 11.723375 170.37206]*10^-5;

% Height correction
a_ht = 2.53*10^-5;
b_ht = 5.49*10^-3;
c_ht = 1.14*10^-3;

% Wet Mapping Function
% See Table5.5
lat_wmf = lat_hmf;
abc_w2po = [5.8021897*10^-4 1.4275268*10^-3 4.3472961*10^-2
5.6794847*10^-4 1.5138625*10^-3 4.6729510*10^-2
5.8118019*10^-4 1.4572752*10^-3 4.3908931*10^-2
5.9727542*10^-4 1.5007428*10^-3 4.4626982*10^-2
6.1641693*10^-4 1.7599082*10^-3 5.4736038*10^-2];

hs_km = h/1000;
DOY_atm = DOY-28;
if lat < doy_atm =" DOY_atm+365.25/2;" lat =" abs(lat);" doyr_atm =" DOY_atm*day2rad;" cost =" cos(DOYr_atm);" a =" abc_avg(1,1);" b =" abc_avg(1,2);" c =" abc_avg(1,3);" i =" 1:4">= lat_hmf(i,1)) & (lat <= lat_hmf(i+1,1)) dlat = (lat-lat_hmf(i))/(lat_hmf(i+1)-lat_hmf(i)); daavg = abc_avg(i+1,1)-abc_avg(i,1); daamp = abc_amp(i+1,1)-abc_amp(i,1); aavg = abc_avg(i,1)+dlat*daavg; aamp = abc_amp(i,1)+dlat*daamp; a = aavg+aamp*cost; % + or - i don't sure dbavg = abc_avg(i+1,2)-abc_avg(i,2); dbamp = abc_amp(i+1,2)-abc_amp(i,2); bavg = abc_avg(i,2)+dlat*dbavg; bamp = abc_amp(i,2)+dlat*dbamp; b = bavg+bamp*cost; % + or - i don't sure dcavg = abc_avg(i+1,3)-abc_avg(i,3); dcamp = abc_amp(i+1,3)-abc_amp(i,3); cavg = abc_avg(i,3)+dlat*dcavg; camp = abc_amp(i,3)+dlat*dcamp; c = aavg+aamp*cost; % + or - i don't sure end end if lat >= lat_hmf(5,1)
a = abc_avg(5,1)+abc_amp(5,1)*cost; % + or - i don't sure
b = abc_avg(5,2)+abc_amp(5,2)*cost; % + or - i don't sure
c = abc_avg(5,3)+abc_amp(5,3)*cost; % + or - i don't sure
end
sine = sin(elev*deg2rad);
cose = cos(elev*deg2rad);

beta = b/(sine+c);
gamma = a/(sine+beta);
topcon = 1+(a/(1+b/(1+c)));
hmf(1,1) = topcon/(sine+gamma);
hmf(2,1) = -topcon*cose/((sine+gamma)^2*...
(1-a/((sine+beta)^2*(1-b/(sine*c)^2))));

beta = b_ht/(sine+c_ht);
gamma = a_ht/(sine+beta);
topcon = 1+a_ht/(1+b_ht/(1+c_ht));
ht_corr_coef = 1/sine-topcon/(sine+gamma);
ht_corr = ht_corr_coef*hs_km;
hmf(1,1) = hmf(1,1)+ht_corr;
dhcc_del = -cose/sine^2+topcon*cose/((sine+gamma)^2*...
(1-a_ht/((sine+beta)^2*(1-b_ht/(sine*c_ht)^2))));
dht_corr_del = dhcc_del*hs_km;
hmf(2,1) = hmf(2,1)+dht_corr_del;
dht_corr_del = dhcc_del*hs_km;
hmf(2,1) = hmf(2,1)+dht_corr_del;
if lat <= lat_wmf(1,1) alat = abc_w2po(1,1); blat = abc_w2po(2,1); clat = abc_w2po(3,1); end for i = 1:4 if (lat >= lat_wmf(i,1)) & (lat <= lat_wmf(i+1,1)) dll = (lat-lat_wmf(i,1))/(lat_wmf(i+1,1)-lat_wmf(i,1)); da = abc_w2po(i+1,1)-abc_w2po(i,1); alat = abc_w2po(i,1)+dll*da; db = abc_w2po(i+1,2)-abc_w2po(i,2); blat = abc_w2po(i,2)+dll*db; dc = abc_w2po(i+1,3)-abc_w2po(i,3); clat = abc_w2po(i,3)+dll*dc; end end if lat >= lat_wmf(5,1)
alat = abc_w2po(5,1);
blat = abc_w2po(5,2);
clat = abc_w2po(5,3);
end
sinelat = sin(elev*deg2rad);
coselat = cos(elev*deg2rad);
betalat = blat/(sinelat+clat);
gammalat = alat/(sinelat+betalat);
topconlat = 1+alat/(1+blat/(1+clat));
wmf(1,1) = topconlat/(sinelat+gammalat);
wmf(2,1) = -topconlat/((sinelat+gammalat)^2*(coselat-alat/...
((sinelat+betalat)^2*coselat*(1-blat/(sinelat*clat)^2))));
hmf(3,1) = wmf(1,1);
hmf(4,1) = wmf(2,1);

hmf;

%%%%%%%%%%END%%%%%%%%%%


.

Lagrange: Satellite Orbit and Satellite Clock Interpolation Source Code for GPS/GPS PPP/GNSS PPP

.

function [xsp3,ysp3,zsp3,clksat,P1,P2,sp3Vx,sp3Vy,sp3Vz ...
,RV1,RV2,RVIO,saas,Elev,AZIMUTH,PH1RV,PH2RV,PHRVIO] ...
= LG_SP3(DOY,rT,PH1,PH2,PR1,PR2,SP3time,SP3X,SP3Y,SP3Z,SATclk,x,SunX,SunY,SunZ,PTRH)

% Lagrange:satellite orbit and satellite clock interpolation
% called by EKF. Computes PHI (dynamic model) for receiver clock
%
% %%%%%%%%%% HELP %%%%%%%%%%
%
% [xsp3,ysp3,zsp3,clksat,P1,P2,sp3Vx,sp3Vy,sp3Vz ...
% ,RV1,RV2,RVIO,saas,Elev,AZIMUTH,PH1RV,PH2RV,PHRVIO] ...
% = LG_SP3(DOY,rT,PH1,PH2,PR1,PR2,SP3time,SP3X,SP3Y,SP3Z,SATclk,x,SunX,SunY,SunZ,PTRH)
%
% Input data
% DOY = day of year
% rT = nominal (RINEX) time for epoch ep [GPSweek GPSsec], [numep x 2]
% PH1 = Carrier phase L1,cycles [numep+1 x numsat]
% PH2 = Carrier phase L2,cycles [numep+1 x numsat]
% PR1 = Pseudorange C/A code, metres [numep+1 x numsat]
% PR2 = Pseudorange P code , metres [numep+1 x numsat]
% SP3time = GPS time of SP3 file [GPSweek GPSsec], [numep x 2]
% SP3X, SP3Y, SP3Z = Satellite coordinate (WGS84), metres [numep+1 x numsat]
% x = Aproximate station coordinates, metres [3 x 1]
% SATclk = Precise clock correction, microseconds [numep+1 x numsat]
% SunX, SunY, SunZ = Sun coordinate (WGS84), metres [numep+1 x numsat]
% PTRH = [P TC RH]
%
% OUTPUT
% xsp3
% ysp3
% zsp3
% clksat
% P1
% P2
% sp3Vx
% sp3Vy
% sp3Vz
% RV1
% RV2
% RVIO
% saas
% Elev
% AZIMUTH
% PH1RV
% PH2RV
% PHRVIO
%
% Modified by Phakphong Homniam
% October 13, 2002
% Original Mathcad source code by Boonsap Witchayangkoon, 2000

%%%%%%%%%% CONSTANT DECLARATION %%%%%%%%%%

% Global constants of the relevant parameters in constant_global.m file
constant_global

global maskangle antenna_offset IIR

% Satellite mask angle
maskangle = 15*deg;

% Satellite Antenna Offset
antenna_offset = [0.279 0 1.023]';

% GPS Block IIR satellite
IIR = 13;
% For GPS Block IIR satellite, IGS has agreed upon to
% the satellite phase center offset [0 0 0]

% Langrange InterPolation Constants
% Number of points used in LG orbit interpolation
nXYZ = 9;
% Number of points used in LG clock interpolation
nCLK = 5;

% Environmental data for this station
P = PTRH(1); % Pressure P in milibar (mb)
TC = PTRH(2); % Temperature TC in degree CELCIUS
RH = PTRH(3); % Humedity H(%)
T = TC+273.2; % Temperature T in Kelvin

%%%%%%%%%% BEGIN %%%%%%%%%%

PRN = PH1(1,:);
PH1(1,:) = []; SP3X(1,:) = [];
PH2(1,:) = []; SP3Y(1,:) = [];
PR1(1,:) = []; SP3Z(1,:) = [];
PR2(1,:) = []; SATclk(1,:) = [];
rT = rT(:,2); SP3time = SP3time(:,2);

x = x(:);
llh = ecef2lla(x',1);
wetzenith = SWZD(P,T,RH);
dryzenith = SHZD(P,T,RH);
[tep,NumSat] = size(PR1);
for j = 1:NumSat
fprintf('\nProcessing satellite number %s\n',num2str(PRN(j)));
CLKsat(1:tep,1) = NAJ;
sp3x(1:tep,1) = NAJ;
sp3y(1:tep,1) = NAJ;
sp3z(1:tep,1) = NAJ;
Prange1(1:tep,1) = NAJ;
Prange2(1:tep,1) = NAJ;
Vx(1:tep,1) = NAJ;
Vy(1:tep,1) = NAJ;
Vz(1:tep,1) = NAJ;
R1(1:tep,1) = NA;
R2(1:tep,1) = NA;
RIO(1:tep,1) = NA;
Tropo(1:tep,1) = NA;
ELEV(1:tep,1) = NA;
azimuth(1:tep,1) = NA;
PH1R(1:tep,1) = NA;
PH2R(1:tep,1) = NA;
PHRIO(1:tep,1) = NA;

r = size(IIR,1);
% Assume no satellite offset
NoOffset = 1;
% NoOffset = 0;
% for h = 1:r
% if PRN(j) == IIR(h)
% NoOffset = 1;
% end
% end
for i = 1:tep
% i
rT(i);
if (PR1(i,j) ~= NA) & (PR2(i,j) ~= NA) & (PH1(i,j) ~= NA) & (PH2(i,j) ~= NA)
CLKsat(i) = SVCLK(rT(i),SP3time,SATclk(:,j) ...
,PR1(i,j),nCLK,NAJ);
if CLKsat(i) ~= NAJ
[X,V] = SatPo(rT(i),PR1(i,j),CLKsat(i) ...
,SP3time,SP3X(:,j),SP3Y(:,j),SP3Z(:,j) ...
,x,nXYZ,[SunX(i),SunY(i),SunZ(i)]',NoOffset,NAJ);
rho = sqrt((X(1)-x(1))^2+(X(2)-x(2))^2+(X(3)-x(3))^2);
elev = elevation(X,x);
if elev < maskangle
continue
end
azimuth(i) = Az(llh(1),llh(2),X-x);
ELEV(i) = elev;
sp3x(i) = X(1);
sp3y(i) = X(2);
sp3z(i) = X(3);
Vx(i) = V(1);
Vy(i) = V(2);
Vz(i) = V(3);
nmf = NeillMF(DOY,llh(1)/deg,llh(3),elev/deg);
Tropo(i) = dryzenith*nmf(1)+wetzenith*nmf(2);
Relativity = -2*dot(X(:),V(:))/c;
h = rho-c*CLKsat(i)*10^-6+Tropo(i)-Relativity;
R1(i) = PR1(i,j)-h;
R2(i) = PR2(i,j)-h;
% See eq.9.46
RIO(i) = cf1*PR1(i,j)-cf2*PR2(i,j)-h;
PH1R(i) = lamda1*PH1(i,j)-h;
PH2R(i) = lamda2*PH2(i,j)-h;
% See eq.9.51
PHRIO(i) = lamda1*(cf1*PH1(i,j)-cf1f2*PH2(i,j))-h;
% if (abs(R1(i)) < rclock) & (abs(R2(i)) < rclock)
Prange1(i) = PR1(i,j);
Prange2(i) = PR2(i,j);
% end
end
end
end
xsp3(1:tep,j) = sp3x(:);
ysp3(1:tep,j) = sp3y(:);
zsp3(1:tep,j) = sp3z(:);
clksat(1:tep,j) = CLKsat(:);
P1(1:tep,j) = Prange1(:);
P2(1:tep,j) = Prange2(:);
sp3Vx(1:tep,j) = Vx(:);
sp3Vy(1:tep,j) = Vy(:);
sp3Vz(1:tep,j) = Vz(:);
RV1(1:tep,j) = R1(:);
RV2(1:tep,j) = R2(:);
RVIO(1:tep,j) = RIO(:);
saas(1:tep,j) = Tropo(:);
Elev(1:tep,j) = ELEV(:);
AZIMUTH(1:tep,j) = azimuth(:);
PH1RV(1:tep,j) = PH1R(:);
PH2RV(1:tep,j) = PH2R(:);
PHRVIO(1:tep,j) = PHRIO(:);
end

%%%%%%%% RESULT %%%%%%%%

xsp3 = [PRN;xsp3];
ysp3 = [PRN;ysp3];
zsp3 = [PRN;zsp3];
clksat = [PRN;clksat];
P1 = [PRN;P1];
P2 = [PRN;P2];
sp3Vx = [PRN;sp3Vx];
sp3Vy = [PRN;sp3Vy];
sp3Vz = [PRN;sp3Vz];
RV1 = [PRN;RV1];
RV2 = [PRN;RV2];
RVIO = [PRN;RVIO];
saas = [PRN;saas];
Elev = [PRN;Elev];
AZIMUTH = [PRN;AZIMUTH];
PH1RV = [PRN;PH1RV];
PH2RV = [PRN;PH2RV];
PHRVIO = [PRN;PHRVIO];

%%%%%%%%%%END%%%%%%%%%%