%Data Sets for pure_data.m-Columns
% 1 point name
% 2 year
% 3 day of year
% 4 hour California time
% 5 minute
% 6 minutes past Jan 8, 2005 00:08:00 UT
% 7 gravity dial reading
% 8 height correction (cm)
% 9 latitude degrees North
% 10 longitude degrees East
% 11 absolute height elliptical (m)
%%%%%%%%%%%%%

load fc05.dat;
load fc04.dat;
load fc05.dat;
load fc04.dat;
name=fc05(:,1);
minutes=fc05(:,6);
dial=fc05(:,7);
hcorr=[fc04(:,8);fc05(:,8)];
lat=[fc04(:,9);fc05(:,9)];
long=[360-fc04(:,10);fc05(:,10)];
height=[fc04(:,11);fc05(:,11)];

cname=fc04(:,1);
cminutes=fc04(:,6);
cdial=fc04(:,7);

%Gravity Conversion
%formula derived elsewhere
%%%%%%%%%%%%%%%%%%%
g_dial=dial.*1.0576-6.6959;
cg_dial=cdial.*1.0576-6.6959;
ga_dial=[cg_dial;g_dial];


%Tidal Correction
%%%%%%%%%%%%%%%%%%
load tides.dat
Atides=[(tides(:,3)-8).*1440+60.*(-8+tides(:,4))+tides(:,5) tides(:,7)./1000];
plot(Atides(:,1),Atides(:,2))
R=floor((minutes)./5+.5)+1;
g_tide=g_dial-Atides(R,2);

load ctides.dat
Btides=[(ctides(:,3)-11).*1440+60.*(-8+ctides(:,4))+ctides(:,5) ctides(:,7)./1000];
figure;
plot(Btides(:,1),Btides(:,2))
S=floor((cminutes)./5+.5)+1;
cg_tide=cg_dial-Btides(S,2);
aname=[cname;name];
ga_tide=[cg_tide;g_tide];

%Drift Correction
%%%%%%%%%%%%%%%%%
%extract Base Camp Values from Tidal correction
%fit linear function
%correct ga_tide according to linear function

%Latitude Correction
%%%%%%%%%%%%%%%%%%%%
lat_corr=978032*(1 + 5.2789e-3*(sin(lat*pi/180)).^2 - 2.35e-6*(sin(lat*pi/180)).^4)-978032*(1 + 5.2789e-3*(sin(max(lat)*pi/180)).^2 - 2.35e-6*(sin(max(lat)*pi/180)).^4);
g_lat=[ga_tide-lat_corr];

%Free Air Correction
%%%%%%%%%%%%%%%%%%%%
aheight=height-hcorr./100;
g_freeair=[g_lat+.307.*(height-min(height))]

%Bouguer Correction
%%%%%%%%%%%%%%%%%%%
g_boug=[g_lat+(.307-.112).*(height-min(height))];

%Terrain Correction
%%%%%%%%%%%%%%%%%%%
%get function

%GRAPHING
%%%%%%%%%
%[x,y]=meshgrid(min(long):.001:max(long), min(lat):.001:max(lat));
[x,y]=meshgrid(245.435:.0001:245.470, 33.99:.0001:34.04);
z=griddata(long,lat,g_boug,x,y);
figure;
shading interp
pcolor(x,y,z);
title('Bouguer gravity');
%[x,y]=meshgrid(min(long):.0001:max(long), min(lat):.0001:max(lat));
%[x,y]=meshgrid(245.435:.0001:245.470, 33.99:.0001:34.04);
shading interp
colorbar;
hold on
plot(long, lat, 'ko')

z=griddata(long,lat,g_freeair,x,y);
figure;
pcolor(x,y,z);
title('Free air gravity');
shading interp
colorbar;
hold on
plot(long, lat, 'ko')

z=griddata(long,lat,ga_tide,x,y);
figure;
pcolor(x,y,z);
title('Tidal corrected gravity');
shading interp
colorbar;
hold on
plot(long, lat, 'ko')

z=griddata(long,lat,height,x,y);
figure;
pcolor(x,y,z);
title('Height');
shading interp
colorbar;
hold on
plot(long, lat, 'ko')