function [NCN,fullind,zeroind]=CNGridSensitivity(Param,source,tolerance,xy,noise)
% Makes a graph of parameters combinations for a spherical source that will
% and will not be detectable by the user-specified gravity network
% inputs:Param = structure with three elements named "Depths","Volumes",
% and "Rhos".
% Param.Depths = a column vector with each of the source depths to
% be investigated
% Param.Volumes = a column vector with each of the source
% volumes to be investigated
% Param.Rhos = a column vector with each of the source
% density contrast to be investigated.
% source=UTM coordinates for source
% tolerance=minimum number of stations for source to be considered
% "detected"
% xy=list of station locations in UTM coords, column vector
% noise=noise level of instrument/survey in mGal
% output:NCN=a matrix containing all combinations used
% fullind=indices for combinations in NCN that are detectable
% zeroind=indices for combinations in NCN that are not detectable
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Code written by Patricia MacQueen (pmacquee@gmail.com). If you use this code, please credit me appropriately!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%variables%%%
Depths=Param.Depths;%depths
Volumes=Param.Volumes;%volumes
rhos=Param.Rhos; %density contrasts
G=6.67e-11;%gravitational constant
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%Cerro Negro (CN) Source%%%%%%%%%%%
%Preallocate N
NCN=zeros(length(Depths)*length(Volumes)*length(rhos),3);
counterind=0;%counter variable to keep track of rows.
goodcountind=0;%holds row locations of detectable parameter combinations
badcountind=0;%holds row locations of undetectable parameter combinations
fullind=zeros(length(Depths)*length(Volumes)*length(rhos),1);%Matrix that will hold detectable combinations of parameters
zeroind=zeros(length(Depths)*length(Volumes)*length(rhos),1);%Matrix that will hold undetectable combinations of parameters
%If you REALLY want to do it in lat long, be my guest....
% %convert lat/lon to local coordinate system with CN as 0,0
% xy(1:length(ggrid),1:2)=0;
% for i=1:length(ggrid)
% xy(i,1:2)=llh2local([ggrid(i,1); ggrid(i,2)],[-86.703339; 12.506016]);
% end
% xy=xy*1000;%convert km to m
%calculate distance to source from each station
r(1:length(xy))=0;
for i=1:length(xy)
r(i)=sqrt((xy(i,1)-source(1))^2+(xy(i,2)-source(2))^2);
end
%Find which stations are above noise level for all depth/volume/rho
%combinations
for i=1:length(Depths)
for j=1:length(Volumes)
for k=1:length(rhos)
delgr=((G*rhos(k)*Depths(i)*Volumes(j))./(r.^2+Depths(i)^2).^(1.5)).*1e5; %calculate gravity at each station in mGal
goodstat=delgr(delgr>noise);%creat a vector containing only gravity values above noise
counterind=counterind+1;
NCN(counterind,1:3)=[Depths(i) Volumes(j) rhos(k)];
if length(goodstat)>=tolerance %if there are more than (tolerance) stations above noise then keep those values, if not, put zeros
goodcountind=goodcountind+1;
fullind(goodcountind)=counterind; %INDEXING IS YOUR FRIEND
elseif length(goodstat)