function [sensitivity]=SensitivityMap(NCN,tolerance,xy,ixy,sxy,noise)
%creates a variable that records how sensitive the user's planned gravity grid is in different
%places in comparison to an "ideal" grid of uniformly spaced stations.
%This code produces an xyz matrix of source locations and corresponding
%sensitivities that can be gridded and displayed as a countour plot in
%Matlab or Surfer.
%inputs: NCN=a gigantic vector generated by GridSensitivity of all the
% possible combinations of depth, volume, and density contrast
% parameters
% tolerance = how many stations that need to have a gravity signal above the noise level
% to consider something "detectable"
% ggrid = UTM coordinates for real station locations ([Northing
% Easting])
% idealgrid = an ideal grid (UTM) of uniformly and densley spaced stations
% that covers the whole survey area with an ideal distribution of stations
% sgrid = a grid of source locations to investigate (UTM)
% noise = noise level below which data is not considered reliable
% (in mGal)
%output: sensitivity= Easting, Northing, sensitivity of grid at that point. This
%value is a comparison of how often "real grid" detects sources at this point,
%and how often the "ideal grid" detects them.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Code written by Patricia MacQueen (pmacquee@gmail.com). If you use this code, please credit me appropriately!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%If you seriously HAVE to do this in lat long, be my guest....
% %convert lat/lon to local coordinate system with center of grid as 0,0
% xy(1:length(ggrid),1:2)=0;
% gridcenter=[mean(ggrid(:,2)) mean(ggrid(:,1))];
% for i=1:length(ggrid)
% xy(i,1:2)=llh2local([ggrid(i,2); ggrid(i,1)],gridcenter);
% end
% xy=xy*1000;%convert km to m
%
% %convert ideal lat/lon to local coordinate system with center of grid as 0,0
% ixy(1:length(idealgrid),1:2)=0;
% for i=1:length(idealgrid)
% ixy(i,1:2)=llh2local([idealgrid(i,2); idealgrid(i,1)],gridcenter);
% end
% ixy=ixy*1000;%convert km to m
%
% %convert source grid lat/lon to local coordinate system with center of grid as 0,0
% sxy(1:length(sgrid),1:2)=0;
% for i=1:length(sgrid)
% sxy(i,1:2)=llh2local([sgrid(i,2); sgrid(i,1)],gridcenter);
% end
% sxy=sxy*1000;%convert km to m
tic%times how long the code takes
realsens=0;%starting counter value for the number of stations detected by the real grid
idealsens=0;%starting counter value for the number of stations detected by the ideal grid
G=6.67e-11;%gravitational constant
sensitivity=zeros(length(sgrid),3);%preallocate matrix for sensitivity values
Fred=waitbar(0,'Ready to WAIT???');%generates a waitbar to let the user know how much longer they have to train for the world championship in thumb-twiddling.
for i=1:length(sxy)
source=sxy(i,:);
%distances to "source"
%real grid
r(1:length(xy))=0;
for j=1:length(xy)
r(j)=sqrt((xy(j,1)-source(1))^2+(xy(j,2)-source(2))^2);
end
clear j
%ideal grid
ir(1:length(ixy))=0;
for j=1:length(ixy)
ir(j)=sqrt((ixy(j,1)-source(1))^2+(ixy(j,2)-source(2))^2);
end
clear j
%real gravity calculations
for m=1:length(NCN)
%real grid gravity
%calculate the gravity signal at the ith source using the mth set
%of parameters and the real grid.
delgr=((G*NCN(m,3)*NCN(m,1)*NCN(m,2))./(r.^2+NCN(m,1)^2).^(1.5)).*1e5;
%ideal grid gravity
%calculate the gravity signal at the ith source using the mth set
%of parameters and the ideal grid.
idelgr=((G*NCN(m,3)*NCN(m,1)*NCN(m,2))./(ir.^2+NCN(m,1)^2).^(1.5)).*1e5;
%if the number of stations in the real grid with a gravity signal above the noise
%level is above tolerance, the real grid at that source location
%gets a "point"
if length(delgr(delgr>noise))>=tolerance
realsens=realsens+1;
end
%if the number of stations in the ideal grid with a gravity signal above the noise
%level is above tolerance, the ideal grid at that source location
%gets a "point"
if length(idelgr(idelgr>noise))>=tolerance
idealsens=idealsens+1;
end
end
%assign a sensitivity value for that source locaiton in the sensitivity
%matrix
sensitivity(i,1:3)=[sgrid(i,1) sgrid(i,2) (realsens/idealsens)*100];
%update the wait bar
waitbar(i/length(sgrid),Fred,'You have to wait THIS much more!! XD');
%reset the counter values for realsens and idealsens to prepare for the
%next source location
realsens=0;
idealsens=0;
end
close(Fred)%closes the waitbar.
toc%stops timing the code and outputs a time value to the screen