%  An example ATG gravity data set
%Elevation(m)  ?g(mgal)   uncertainty(mgal)    Floor
data=[
     0          0            0.02               0
    4.6       -1.00          0.02               1
    8.5       -1.88          0.03               2
   12.4       -3.01          0.03               3
   16.3       -4.13          0.04               4
   20.5       -5.41          0.04               5
   24.4       -6.56          0.06               6
   28.3       -7.72          0.06               7];



% for the second part of the lab, here are the data:


A=[-122.333830   47.416332      102.40     4767.80
-122.333340   47.438332      103.60     4767.80
-122.332660   47.459999       69.40     4776.71
-122.332660   47.472668      109.70     4768.54
-122.332500   47.488834      127.40     4762.91
-122.332660   47.503334      141.40     4751.21
-122.332830   47.512165      130.70     4745.30
-122.323830   47.521332       11.20     4762.80
-122.321670   47.528332        4.20     4756.60
-122.321330   47.536335        3.60     4750.10
-122.328330   47.544834        4.20     4744.40
-122.332830   47.557667        4.50     4735.60
-122.333000   47.575832        4.20     4730.85
-122.332660   47.596001        4.20     4727.96
-122.341670   47.612835       45.10     4719.30
-122.333340   47.625832        8.20     4727.80
-122.325000   47.639999       31.30     4723.85
-122.326670   47.655834       36.20     4725.56
-122.327160   47.672668       84.40     4718.75
-122.335170   47.691002       64.90     4729.24
-122.343330   47.705166      100.20     4727.32
-122.343500   47.723499      132.20     4729.16
-122.343830   47.734665      146.60     4730.46
-122.344330   47.752335      142.00     4741.21
-122.344670   47.763332      145.90     4746.90
-122.344670   47.777832      124.30     4760.00
-122.334170   47.794666      104.80     4768.71];


lon=A(:,1);
lat=A(:,2);
elev=A(:,3);
raw_grav=A(:,4);

% remember to add 976000 to the data to have absolute gravity in mgal

% I gave you latitudes and longitudes, you need to convert to distances in
% km north of thestart of the profile.  You can use the m_map function
% previously used.
%  This returns the distances between the adjacent points using the nominal
%  longitude.
dx=m_lldist(-122.34*ones(size(lat)),lat);
%to get the distance from the starting point (at the south end)
% cumsum adds the dx's sequentially. 
dist=[0 ;cumsum(dx)];

% note that the raw_grav is in mgal 
% I provide a function to give you the gravity predicted by the
% International Gravity Formul (but it returns g in gravity units)
%  g=grav_lat(latitude)

% Add code here to calculate the IGF gravity and substract it from the
% observed gravity, then 

% Add a couple lines to calculate the Free Air Anomaly

% Add line here to calculate the Bouguer Anomaly


% To check your corrections, for the first observation:
% -122.333830   47.416332      102.40     4767.80
% the Free Air anomaly is:     -381.6 gu
% and the Bougier anomaly is:  -494.5 gu
% These are in gu units
% you need to divide these by 10 to get anomalies in mgal

% Add lines here to plot the various anomalies vs location

% Also make a map showing where the data are taken.

%Since the m_map data base has a course coastline, I am providing a higher
%resolution data base of coastlines for the Cascadia area and a different
%function to plot the features.  The data base is in a file called
%Puget_big.mat - it is on the server.  The following commands will plot
%the outline (you need to adjust the longitude and latitude ranges to
%optimize the appearance of your maps.  Note that this makes use of a
%MATLAB feature calle data structures (like objects in C++). args is a data
%structure that contains at least two items (the name of the file
%containing the data and the map limits).

args.file='puget_big';
lon_range=[-125 -120];  %adjust this for better visualization
lat_range=[46 49];         %adjust this for better visualization
args.map_limits=[ lon_range lat_range ];
figure(1);
mapp(args);
hold on
plot(lon,lat,'o')
hold off

