%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%gravity.m
%this is the main function for data reduction of the gravity survey.
%its main steps: 
%               1)gain data from the given global data matrix.
%               2) corrections: (a) tides (b)drift (c) latitude 
%                 (d) free-air (e) Bouguer
%               3) plots: anomalies vs. heights; contour map of resulting gravity anomalies 
%
% input: loaded data matrix, and a range (number of last matrix line to use)
%   input file format:
%   Station_number 
%   dial_riading   
%   time(hour) 
%   time (minute)    
%   lat(deg)
%   lat (dec.minutes)
%   longitude (deg)
%   longitude(dec.minutes)
%   elevation (relative to base-camp)
%   GPS antena height (relative to gravimeter)
%   day of measurements (1 being the first, and then consecutively)
%   example:
%   001	  3052.71 10 56 34 2.2493 114 32.591 81.433 1.44 1
%
%   if u change the input file format, u only need to change the first
%   section in this file (data retriving), and the section relating to the
%   base_camp data in "drift_corr.m"
%
%   output: lonitudes, latitudes, heights, resulting gravity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%function [longs, lats, heights, result] = gravity (data_matrix, range)
range = 188+24;
%range = 29;
load ('grav_all');

data_matrix = grav_all;

%get columns from the data matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

lat_degrees = data_matrix (1:range, 5);
lat_minutes = data_matrix (1:range, 6)./60;
longs = -114 - data_matrix(1:range, 8)./60;
lats = lat_degrees + lat_minutes;
days = data_matrix (1:range, 11);
hours = data_matrix(1:range, 3);
minutes = data_matrix (1:range,4 );
dial_readings = data_matrix(1:range, 2);
heights = data_matrix(1:range, 9) - data_matrix(1:range, 10) + 1.293;


%data reduction steps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%g_after_tidal = tidal_correction (dial_readings, days, hours, minutes);
%dial_after_corr = drift_corr(g_after_tidal, days, hours, minutes, data_matrix((range+1):length(data_matrix(:,1)), :));
%raw_g = dial2gravity (dial_after_corr);
raw_g = dial2gravity (dial_readings);

%raw_g = 1.0579*dial_after_corr';


lat_radians = (lat_degrees + lat_minutes)*pi/180;
g_after_lat = lat_correction (raw_g, lat_radians);
g_after_car= g_after_lat - g_after_lat(1); 
g_after_faa = g_after_car + 0.307.*heights';
g_after_bouguer = bouguer_corr (g_after_faa, heights');
result = g_after_bouguer;

mean_g = mean(result);


%plotting of results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


figure;
plot (heights, result-mean_g, 'ks');
title ('after bouguer');
grid on;
hold on;
for i=1:range
	plot(heights(i), result(i)-mean_g, '.w', 'markersize', 20);
	text( heights(i), result(i)-mean_g, sprintf(' %d',data_matrix(i,1)));
end;
%%%%%%%%
figure;
plot (heights, g_after_faa, 'ks', 'markersize', 14);;
title ('after faa');
grid on;
hold on;
for i=1:range
	plot(heights(i), g_after_faa(i), '.w', 'markersize', 20);
	text( heights(i),g_after_faa(i),sprintf(' %d',data_matrix(i,1)));
end;
%%%%%%%%%5
%figure;
%plot (heights, g_after_tidal-mean(g_after_tidal), 'ks');
%title ('tidal correction');
%grid on;
%hold on;
%for i=1:range
%	plot(heights(i), g_after_lat(i)-mean(g_after_lat)  , '.w', 'markersize', 20);
%	text( heights(i), g_after_lat(i)-mean(g_after_lat) , sprintf(' %d',data_matrix(i,1)));
%end;
%%%%%%%%%%%%

figure;
plot (heights, g_after_lat-mean(g_after_lat), 'ks');
title ('lat corr');
grid on;
hold on;
for i=1:range
	plot(heights(i), g_after_lat(i)-mean(g_after_lat)  , '.w', 'markersize', 20);
	text( heights(i), g_after_lat(i)-mean(g_after_lat) , sprintf(' %d',data_matrix(i,1)));
end;
%%%%%%%%%%%%


[x,y] = meshgrid (  min(longs) : 0.0001 : max(longs),  min(lat_degrees + lat_minutes) : 0.0001 : max(lat_degrees + lat_minutes));

z = griddata (longs, lat_degrees + lat_minutes, g_after_bouguer', x, y,'cubic');
figure;
%contourf (x,y,z,30);
pcolor (x, y, z);
colorbar;
shading  interp;
title 'g after Bouguer , Lat and FAA corrections');
grid on;
hold on;
for i=1:range
	plot(longs(i), lat_degrees(i)+lat_minutes(i), '.w', 'markersize', 20);
	text(longs(i), lat_degrees(i)+lat_minutes(i), sprintf(' %d',data_matrix(i,1)));
end;
axis tight
%%%%%%%%%%%%

figure;
contourf(x,y,z,30);
title 'g after Bouguer , Lat and FAA corrections');
colorbar;
hold on;
for i=1:range
	plot(longs(i), lat_degrees(i)+lat_minutes(i), '.w', 'markersize', 20);
	text(longs(i), lat_degrees(i)+lat_minutes(i), sprintf(' %d',data_matrix(i,1)));
end;
axis tight
%%%%%%%%%%%%


figure;
plot3(longs, lat_degrees+lat_minutes, heights,'.-k');
hold on;
grid on;

for i=1:range
	plot3(longs(i), lat_degrees(i)+lat_minutes(i), heights(i),'.w', 'markersize', 20);
	text(longs(i), lat_degrees(i)+lat_minutes(i), heights(i), sprintf(' %d',data_matrix(i,1)));
end;

% make input for g_o_d.m
gravity_c(:,1) = longs;
gravity_c(:,2) = lats;
gravity_c(:,3) = heights;
gravity_c(:,4) = g_after_bouguer';
save gravity_c