%free air anomaly GRS80,free air anomalies reduction to EIGEN-CG03C,free air anomalies reduction to model of remaining topography
load grav67.dat; %Point ID,Latitude(deg),Longitude(deg),h(m),DgF(GRS67)(mGal)
load grav80eig.txt; %Point ID,Latitude(deg),Longitude(deg),h(m),DgGM(EIGEN-CG03C)(mGal)
load grav80rtm.txt; %Point ID,Latitude(deg),Longitude(deg),h(m),DgRTM(Residual Terrain Model)(mGal)
format long
GRS80=grav67;
EIGEN=grav67;
MRT=grav67;
for i=1:1173;
g0GRS67(i,1)=(9.78031846*(1+(0.0053024*(sin(grav67(i,2)*pi/180))^2)-(0.0000059*(sin(2*grav67(i,2)*pi/180))^2)))*(10^5); %Normal Gravity GRS67 (mGal)
g0GRS80(i,1)=(9.78032677*((1+(0.001931851353*sin(grav67(i,2)*pi/180)^2))/((1-(0.0066943800229*sin(grav67(i,2)*pi/180)^2))^(1/2))))*(10^5); %Normal Gravity GRS80 (mGal)
GRS80(i,5)=grav67(i,5)+g0GRS67(i,1)-g0GRS80(i,1); %Calculation of DgF(GRS80) (mGal)
EIGEN(i,5)=GRS80(i,5)-grav80eig(i,5); %Calculation of free air anomalies reduction to EIGEN-CG03C
MRT(i,5)=EIGEN(i,5)-grav80rtm(i,5); %Calculation of free air anomalies reduction to model of remaining topography
DgFGRS80(i,1)=GRS80(i,5); %Free air anomalies DgF(GRS80) table
DgEIGEN(i,1)=EIGEN(i,5); %Free air anomalies EIGEN-CGO3C table
DgRES(i,1)=MRT(i,5); %Free air anomalies MRT table
end
fid=fopen('GRS80.dat','w');
for i=1:1:1173
for j=1:1:5
fprintf(fid,'%15.3f',GRS80(i,j));
end
fprintf(fid,'\n');
end
fclose(fid);
fid=fopen('GRS80_EIGEN.dat','w');
for i=1:1:1173
for j=1:1:5
fprintf(fid,'%15.3f',EIGEN(i,j));
end
fprintf(fid,'\n');
end
fclose(fid);
fid=fopen('MRT.dat','w');
for i=1:1:1173
for j=1:1:5
fprintf(fid,'%15.3f',MRT(i,j));
end
fprintf(fid,'\n');
end
fclose(fid);
save ANSWERS123 GRS80 EIGEN MRT %Saving the results
%mean value,max,min,standard deviation calculation of three fields
GRS80mean=mean(DgFGRS80);
EIGENmean=mean(DgEIGEN);
RESmean=mean(DgRES);
GRS80max=max(DgFGRS80);
EIGENmax=max(DgEIGEN);
RESmax=max(DgRES);
GRS80min=min(DgFGRS80);
EIGENmin=min(DgEIGEN);
RESmin=min(DgRES);
GRS80std=std(DgFGRS80);
EIGENstd=std(DgEIGEN);
RESstd=std(DgRES);
save ANSWERS4 GRS80mean EIGENmean RESmean GRS80max EIGENmax RESmax GRS80min EIGENmin RESmin GRS80std EIGENstd RESstd %Saving the results
%Calculation of analytic covariance function from empirical covariance function
load grav80eigenrtmcov.txt %(deg),C(DgF,DgF)
for i=1:8;
stat1(i,1)=grav80eigenrtmcov(i,1);
C(i,1)=grav80eigenrtmcov(i,2);
end
n=0;
mo=0;
mi=0;
k=sum(stat1);
l=sum(log(C))/8;
m=0;
for i=1:8;
n=n+(stat1(i,1)*log(C(i,1)));
m=m+(stat1(i,1))^2;
mo=mo+(log(C(i,1)))^2;
mi=mi+(log(C(i,1)));
end
b=((n-(k*l))/(m-(k^2/8))); %Calculation of b
a=exp(l-(b*k/8)); %Calculation of a
for i=1:36;
stat(i,1)=grav80eigenrtmcov(i,1);
Cemp(i,1)=grav80eigenrtmcov(i,2);
Canal(i,1)=a*exp(b*stat(i,1));
Empanal(i,1)=stat(i,1); %Saving empirical and analytic function in a table named "Empanal"
Empanal(i,2)=Cemp(i,1);
Empanal(i,3)=Canal(i,1);
end
r2=((n-(k*l))^2/((m-(k^2/8))*(mo-(mi^2/8)))); %Calculation of r^2
r=sqrt(r2); %Calculation of r
plot(stat,Cemp,'r')
hold on
plot(stat,Canal,'b')
legend ('Empirical', 'Analytic')
title('Empirical and analytic covariance function')
xlabel('Spherical distance (degree)')
ylabel('Covariance function (10 ^-5 ms^-2)^2')
hold off
save covarience Empanal r2 r a b %Saving "Empanal,r2
%Forecast of free air anomaly MRT to P point(=59.405,=-116.333)
fp=(59.405)*pi/180;
lp=(-116.333)*pi/180;
for i=1:1173
adis(i,1)=MRT(i,1);
adis(i,2)=MRT(i,2)*pi/180;
adis(i,3)=MRT(i,3)*pi/180;
adis(i,4)=MRT(i,4);
adis(i,6)=MRT(i,5);
k1=sin(adis(i,2))*sin(fp);
k2=cos(adis(i,2))*cos(fp);
k3=cos(adis(i,3))*cos(lp);
k4=sin(adis(i,3))*sin(lp);
adis(i,5)=acos(k1+(k2*(k3+k4)));
end
% BUBBLE SORT OPTIMIZATION
for i=1:1172
for j=1:1173-i
if (adis(j,5)>adis(j+1,5))
temp=adis(j,6);
adis(j,6)=adis(j+1,6);
adis(j+1,6)=temp;
temp=adis(j,5);
adis(j,5)=adis(j+1,5);
adis(j+1,5)=temp;
temp=adis(j,4);
adis(j,4)=adis(j+1,4);
adis(j+1,4)=temp;
temp=adis(j,3);
adis(j,3)=adis(j+1,3);
adis(j+1,3)=temp;
temp=adis(j,2);
adis(j,2)=adis(j+1,2);
adis(j+1,2)=temp;
temp=adis(j,1);
adis(j,1)=adis(j+1,1);
adis(j+1,1)=temp;
end
end
end
for i=1:5
adist(i,1)=adis(i,1);
adist(i,2)=adis(i,2)*180/pi;
adist(i,3)=adis(i,3)*180/pi;
adist(i,4)=adis(i,4);
adist(i,5)=adis(i,5);
end
for i=1:5
for j=1:5
k1=sin(adist(i,2)*pi/180)*sin(adist(j,2)*pi/180);
k2=cos(adist(i,2)*pi/180)*cos(adist(j,2)*pi/180);
k3=cos(adist(i,3)*pi/180)*cos(adist(j,3)*pi/180);
k4=sin(adist(i,3)*pi/180)*sin(adist(j,3)*pi/180);
nearij(i,j)=acos(k1+(k2*(k3+k4))); %Calculation of nearest i,j points
CDgiDgj(i,j)=a*exp(b*nearij(i,j)*180/pi); %Calculation of C(Dgi,Dgj)
end
end
for i=1:5
Dgires(i,1)=adis(i,6); %Vector of DgMRT
CDgpDgi(i,1)=a*exp(b*adist(i,5)*180/pi); %Vector of C'(Dgp,Dgi)
end
CDgiDgj1=inv(CDgiDgj); %Matrix CDgiDgj inverse
Dgp1=CDgpDgi'*(CDgiDgj1*Dgires); %Forecast of free air anomaly MRT to P point
CDgpDgi=CDgpDgi';
CDgiDgje=CDgiDgj;
for i=1:5
CDgiDgje(i,i)=CDgiDgje(i,i)+2; %Matrix CDgiDgj with errors
end
CDgiDgj2=inv(CDgiDgje); %CDgiDgje inverse
Dgp2=CDgpDgi*(CDgiDgj2*Dgires); %Forecast of free air anomaly MRT to P point with errors
forcerror1=sqrt(a-(CDgpDgi*CDgiDgj1*CDgpDgi'));
forcerror2=sqrt(a-(CDgpDgi*CDgiDgj2*CDgpDgi'));
dg80=Dgp1+16.532-1.854;
dg80e=Dgp2+16.532-1.854;
save forecast adist nearij CDgpDgi CDgiDgj CDgiDgj1 Dgp1 CDgiDgj2 CDgiDgje forcerror1 forcerror2 Dgires Dgp2 dg80 dg80e
%least squares collocation
pcoll=adist;
for i=1:5
pcoll(i,6)=Dgires(i,1);
dist(i,1)=adist(i,5);
end
g=980000;
R=6371000;
A=425;
s=0.999617;
t=cos(nearij);
t1=cos(dist);
for i=1:5
L1(i,1)=sqrt(1-(2*s*t1(i,1))+s^2);
N1(i,1)=1+L1(i,1)-(s*t1(i,1));
CNDG(i,1)=(A/g)*R*s*(log(2/N1(i,1))-(t1(i,1)*s));
for j=1:5
L(i,j)=sqrt(1-(2*s*t(i,j))+s^2);
N(i,j)=1+L(i,j)-(s*t(i,j));
CII(i,j)=A*s^2*((1/L(i,j))-1-(log(2/N(i,j))));
end
end
CII1=inv(CII);
b1=CII1*Dgires;
CNDG=CNDG';
NGEOID=CNDG*b1;
geoidfin=NGEOID-19.608+0.326;
save collocation NGEOID CNDG b1 CII CII1 geoidfin
