ClickSim.asv

function [RLforHist,pDetTotal] = ...
    ClickSim(sortedTLVec,diveDepth_mean,diveDepth_std,SL_mean,SL_std,...
    descAngle_mean,descAngle_std,clickStart_mean,clickStart_std,...
    directivity,minAmpSide_mean,...
    minAmpBack_mean,botAngle_std,descentPerc)
%Based on Kait Frasier ClickMethod
%
RLforHist = [];
for itr_n = 1:n  % number of simulations loop
    if rem(itr_n,100) == 0
        fprintf('TL computation %d of %d\n', itr_n, n)
    end
    %%%%% Location Computation %%%%%
    % rand location
    randVec = ceil(rand(2,N)'.*repmat([2*maxRange, 2*maxRange], [N, 1]))...
        - repmat([maxRange, maxRange], [N, 1]);
    [theta, rho] = cart2pol(randVec(:,1),randVec(:,2));  % convert to polar coord.
    clear randVec  % trying to save on memory
    
    % trim out the locations that are beyond the max range (corners of the
    % 2*maxRange X 2*maxRange square, since now we are using a pi*maxRange^2
    % circle)
    jjj = 1;
    rho2 = [];
    theta2 = [];
    for iii = 1:length(rho)
        if rho(iii) < maxRange
            rho2(jjj,1) = rho(iii);
            theta2(jjj,1) = theta(iii);
            jjj = jjj+1;
        end
    end
    thetaDeg = 180 + (theta2*180/pi);
    clear theta rho
    
    % go from angle to ref indices
    [angleRef,radRef] = angle_ref_comp(thetaDeg,rho2,thisAngle);
    
    %%%%% Depth Computation %%%%%
    % Compute bottom depth at each randomly selected point
    count0 = 1;
    tempDepth = zeros(size(angleRef));
    keepPoint = ones(size(angleRef));
    
    diveDepthRef = diveDepth_mean(itr_n) + diveDepth_std(itr_n)...
        *randn(size(angleRef)); % add variation to dive depth
    
    % If there are whales below the seafloor, place them above
    burrowingWhaleIdx = find(diveDepthRef>=sd);
    while ~isempty(burrowingWhaleIdx)
        diveDepthRef(burrowingWhaleIdx) = diveDepth_mean(itr_n)...
            + diveDepth_std(itr_n)*randn(size(burrowingWhaleIdx)); % add variation to dive depth,
        burrowingWhaleIdx = find(diveDepthRef>=sd);
    end
    % Remove unwanted points from the body of points that will be run
    % through the rest of the model.
    rho2 = rho2(keepPoint == 1);
    theta2 = theta2(keepPoint == 1);
    thetaDeg = thetaDeg(keepPoint == 1);
    radRef = radRef(keepPoint == 1);
    angleRef = angleRef(keepPoint == 1);
    
    % Assign last n% to a descent phase
    % Choose a depth between start of clicking and destination depth
    % determine off-axis angle
    descentIdx = (floor((1-descentPerc(itr_n,1))*length(rho2))+1:length(rho2))';
    dFactor = rand(size(descentIdx));
    clickStartVec = clickStart_mean(itr_n,1) + clickStart_std(itr_n,1).*randn(size(descentIdx));
    
    % If there are whales above the sea surface, put them below it.
    flyingWhaleIdx = find(clickStartVec<1);
    while ~isempty(flyingWhaleIdx)
        clickStartVec(flyingWhaleIdx) = clickStart_mean(itr_n,1) + clickStart_std(itr_n,1).*randn(size(flyingWhaleIdx));
        flyingWhaleIdx = find(clickStartVec<1);
    end
    descentDelta = dFactor.* (diveDepthRef(descentIdx,:) - clickStartVec);
    diveDepthRef(descentIdx,1) = clickStartVec + descentDelta;
    
    %%%%% Beam Angle Computation %%%%%
    % Assign random beam orientation in horizontal (all orientations equally likely)
    randAngleVec = ceil(rand(size(rho2)).*359);
    % Compute vertical component of shift between animal and sensor (sd =
    % sensor depth)
    dZ = abs(sd - diveDepthRef);
    
    zAngle_180 = ceil(abs(atand(dZ./radRef))+ (botAngle_std(itr_n,1)*randn(size(dZ))));
    % assign descent angle to descending portion
    zAngle_180(descentIdx,1) = ceil(abs(atand(dZ(descentIdx,:)./radRef(descentIdx,:))) -...
        descAngle_mean(itr_n,1) + (descAngle_std(itr_n,1).*randn(size(descentIdx))));
    
    zAngle = make360(zAngle_180); % wrap
    % clear zAngle_180
    
    %%%%% Transmission loss (TL) Computation %%%%%
    % Note, due to computation limitations, directivity does not vary by individual.
    % The beam pattern is considered to be the same for all individuals within an iteration.
    % Compute beam pattern:
    [beam3D,~] = odont_beam_3D(directivity(itr_n,1), [minAmpSide_mean(itr_n,1),minAmpBack_mean(itr_n,1)]);
    
    % Compute variation to add to source level
    SL_adj = SL_std(itr_n,1)*randn(size(zAngle));
    
    RL = nan(size(thetaDeg));
    isheard = zeros(size(thetaDeg));
    %%%%% Transmission Loss Loop %%%%%
    for itr2 = 1:length(thetaDeg)
        % Using vertical and horizontal off axis components, compute beam
        % related transmission loss
        beamTL = beam3D(zAngle(itr2), randAngleVec(itr2));
        
        % Compute location of this animal in the transmission loss matrix:
        % Find which row you want to look at:
        thisRd = rd_all{angleRef(itr2)};
        [~,thisDepthIdx] = min(abs(thisRd - round(diveDepthRef(itr2))));
        
        % record the distance related portion of this transmission loss
        thisSortedTL = real(sortedTLVec{angleRef(itr2)});
        distTL = thisSortedTL(thisDepthIdx,ceil(radRef(itr2)./rr_int));
        
        % Add up all the sources of TL
        RL(itr2,1) = SL_mean(itr_n,1) + SL_adj(itr2) - beamTL - distTL;
        
        % Is the total TL less than the maximum allowed?
        if RL(itr2,1)>=thresh
            isheard(itr2,1) = 1; % detected it
        end
    end
    
    pDetTotal(itr_n,1) = sum(isheard)./length(isheard)';
    detVsLoc = [thetaDeg, rho2, isheard];
    totalSim = rho2';
    detSim = rho2(isheard==1)';
    RL_keep = RL(isheard==1);
    %                      RLforHist(itr_n,:) = histc(RL_keep,RLbins); % made bins go 0.5 to 1.5
    RLforHist(itr_n,:) = hist(RL_keep,RLbins); % move to integer
    
    % Compute detections in range bins, so you can make a histogram if desired
    % Makes more sense for click-based model
    % preallocate
    binTot = zeros(length(binVec)-1,1);
    binDet = zeros(length(binVec)-1,1);
    
    for itr3 = 1:length(binVec)-1
        binTot(itr3) = length(find(totalSim>binVec(itr3) & totalSim<binVec(itr3 +1)));
        binDet(itr3) = length(find(detSim>binVec(itr3) & detSim<binVec(itr3 +1)));
    end
    thisPercent = binDet./binTot;
    % save the bin counts to the overall set, so you can get means and variances per bin.
    binnedPercDet(itr_n,:) = thisPercent';
    binnedCounts(itr_n,:) = binDet';
    
end
%
end