%% 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:00:00 UT
% 7 gravity dial reading
% 8 height correction (cm)
% 9 latitude degrees North
% 10 longitude degrees East
% 11 absolute height elliptical (m)
%%%%%%%%%%%%%

load fc08.dat;
load fc05.dat;
load fc04.dat;
load fc08.dat;
load fc05.dat;
load fc04.dat;

load g_drift.txt;
load bg_drift.txt;
load cg_drift.txt;

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

bname=fc05(:,1);
bminutes=fc05(:,6);
bdial=fc05(:,7);

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

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


%Tidal Correction
%%%%%%%%%%%%%%%%%%
% removed -8 on 1/24/08 from the tides08(:,4)
% load tides08.dat
% Atides=[(tides08(:,3)-9).*1440+60.*(tides08(:,4))+tides08(:,5) tides08(:,7)/1000];
% %plot(Atides(:,1),Atides(:,2))
% R=floor((minutes)./5+.5)+1;
% g_tide=g_dial-Atides(R,2);
% 
% load tides05.dat
% Btides=[(tides05(:,3)-8).*1440+60.*(tides05(:,4))+tides05(:,5) tides05(:,7)/1000];
% figure;
% %plot(Btides(:,1),Btides(:,2))
% S=floor((bminutes)./5+.5)+1;
% bg_tide=bg_dial-Btides(S,2);
% 
% load tides04.dat
% Ctides=[(tides04(:,3)-11).*1440+60.*(tides04(:,4))+tides04(:,5) tides04(:,7)/1000];
% figure;
% %plot(Ctides(:,1),Ctides(:,2))
% T=floor((cminutes)./5+.5)+1;
% cg_tide=cg_dial-Ctides(T,2);

aname=[cname;bname;name];
ga_drift=[cg_drift;bg_drift;g_drift];

%Drift Correction
%%%%%%%%%%%%%%%%%
%extract Base Camp Values from Tidal correction
%fit linear function
%correct ga_tide according to linear function
%performed in excel manually
%corrected between data sets of different years as well

%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_drift-lat_corr];

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

%Bouguer Correction + Free Air Correction
%%%%%%%%%%%%%%%%%%%
% density of 2.5 g/cc
g_boug=[g_lat+(.307-2.5*0.04193).*(height-min(height))];

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


%%%%%%%%%%%%%%%%%%
g_anom=g_boug-g_boug(215)*ones(302,1) ; %bouger value from base station 2008

%GRAPHING
%%%%%%%%%
%%
%[x,y]=meshgrid(min(long):.001:max(long), min(lat):.001:max(lat));
%[x,y]=meshgrid(245.410:.0001:245.470, 33.99:.0001:34.08);
ll = [245.430 33.998]; Dl = 0.05; dll = [Dl*cos(ll(2)*pi/180) Dl];
Rlng = [ll(1):0.0001:ll(1)+dll(1)]; Rlat=[ll(2):0.0001:ll(2)+dll(2)];
[x,y]=meshgrid(Rlng, Rlat);
z=griddata(long,lat,g_anom,x,y);
%z=griddata(long,lat,height,x,y);

figure(1); hold off; 
shading interp
pcolor(x,y,z);
title('Bouguer Anomaly Relative to Base Camp (Drift Corrected)');
%title('Height');

%[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, 'k+')
daspect([cos(ll(2)*pi/180) 1 1] );
fprintf('Range %10.5f %10.5f %10.5f %10.5f\n',ll(1)-360, ll(1)+dll(1)-360, ll(2),ll(2)+dll(2));
hold off

%%

%GRAPHING
%%%%%%%%%
%[x,y]=meshgrid(min(long):.001:max(long), min(lat):.001:max(lat));
%[x,y]=meshgrid(245.410:.0001:245.470, 33.99:.0001:34.08);
%[x,y]=meshgrid(245.435:.0001:245.468, 33.999:.0001:34.04);
% Use same grid as before
z=griddata(long,lat,g_boug,x,y);
figure(2);
shading interp
pcolor(x,y,z);
title('Total Bouguer Gravity (Drift Corrected)');
%[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')
daspect([cos(ll(2)*pi/180) 1 1] );
caxis([3247 3253.2])
hold off

%% Now do the heights
z=griddata(long,lat,height,x,y);
figure(3);
pcolor(x,y,z);
title('Ellipsoidal Height');
shading interp
colorbar;
hold on
plot(long, lat, 'k+')
daspect([cos(ll(2)*pi/180) 1 1] );
hold off

%% Now do a zoom plot along profile
%  245.4473  245.4588   34.0207   34.0346

llz = [245.447 34.0207]; Dlz = 0.015; dll = [Dlz*cos(llz(2)*pi/180) Dlz];
Rlngz = [llz(1):0.0001:llz(1)+dll(1)]; Rlatz=[llz(2):0.0001:llz(2)+dll(2)];
[xz,yz]=meshgrid(Rlngz, Rlatz);

z=griddata(long,lat,g_anom,xz,yz);
figure(4);
shading interp
pcolor(xz,yz,z);
title('Bouguer Anomaly Along Seismic Line (Drift Corrected)');
%[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, 'k.')
daspect([cos(llz(2)*pi/180) 1 1] );
fprintf('Zoom Range %10.5f %10.5f %10.5f %10.5f\n',llz(1)-360, llz(1)+dll(1)-360, llz(2),llz(2)+dll(2));
hold off

% 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')
%