function [p_not p_notg]=add_bubble_interactions_dist(BubbleData,Signature, ArrayGeo,cw,rhow_inf,dt,r)
 % Ozkan Sertlek, 2015, University of Leiden
% adds the bubble interactions to the notional signatures 

%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.see <http://www.gnu.org/licenses/>.

% convert from structure to single variables..
N=ArrayGeo.N;                               % number of guns
Time_Ax=Signature.Time_Ax(1,:);     % time axis
p=Signature.p_unperturbed;                       % notional signature 
x=ArrayGeo.x;                               % x location of gun
y=ArrayGeo.y;                               % y location of gun
zg=ArrayGeo.zg;                             % z location of gun
zr=ArrayGeo.zr;                             % z location of receiver
Ro=BubbleData.Ro(:,1:end-1);         % bubble radius
Uo=BubbleData.Uo(:,1:end-1);         % bubble wall velocity
Pao=BubbleData.P_bubblewall;         % pressure at bubble wall 

% unit conversions
atm_to_Pa=101325;
P_atm=1*atm_to_Pa;
g=9.80665 ;
n=7;B=3000*atm_to_Pa;

% initialize
p_notg(1:N,:)=0*p(1:N,:);
for j=1:N
    j
    h(1:N, 1:length(Time_Ax))=0; 
    for i=1:N
            dx=x(i)-x(j);
            dy=y(i)-y(j);
            dz=zg(i)-zg(j);
            dz_g=zg(i)+zg(j);
            Di=sqrt(dx^2 + dy^2 +dz^2 ); %distance between gun
            Dgi=sqrt(dx^2 + dy^2 +dz_g^2); %distance between ghos
            td_gun=(Di-1)/cw;
            td_ghost=(Dgi-1)/cw;
            [v id_gun]=min(abs(td_gun*1000-Time_Ax));
            [v id_ghost]=min(abs(td_ghost*1000-Time_Ax));
            ind_ght=id_ghost:length(p(i,:));
            ind_src=id_gun:length(p(i,:));
            P_inf=P_atm +rhow_inf*g*zg(j);
            Hi(i,:)=n*(P_inf+B).*( ((Pao(i,:)+B)./(P_inf+B)).^((n-1)/n) -1)./(rhow_inf*(n-1));
            hs(i,1:id_gun-1)=0;
            hg(i,1:id_ghost-1)=0;
            hs(i,ind_src)=rhow_inf.*Ro(i,ind_src-id_gun+1).*(Hi(i,ind_src-id_gun+1)+0.5*(Uo(i,ind_src-id_gun+1).^2))./Di;  %notional pressure
            hg(i,ind_ght)=rhow_inf.*Ro(i,ind_ght-id_ghost+1).*(Hi(i,ind_ght-id_ghost+1)+0.5*(Uo(i,ind_ght-id_ghost+1).^2))./Dgi;  %notional pressure
            if i==j || Di==0
                h(i,:);
            else
                h(i,:)=h(i,:)+hs(i,:)-hg(i,:);
            end
            P_ieff(j,:)=P_inf+h(i,:);
            Heff(j,:)=n*(P_ieff(j,:)+B).*( ((Pao(j,:)+B)./(P_ieff(j,:)+B)).^((n-1)/n) -1)./(rhow_inf*(n-1));
            p_not(j,:)=rhow_inf.*Ro(j,:).*(Heff(j,:)+(Uo(j,:).^2)./2);  %notional pressure
    end
    % adding ghost
     D1=sqrt((zr(1)-zg(j))^2 + r^2);
     D2=sqrt((zr(1)+zg(j))^2 + r^2);
     R0=sqrt((zr(1))^2 + r^2);
     time_delay= (D2-D1)/cw;              
     [v idg]=min(abs(time_delay*1000-Time_Ax)); % index is shifted 
     ind_ght=idg:(length(p_not(j,:))+idg-1);
     P_image(1:(length(p_not(j,:))+idg))=0;
     P_image(1:idg)=0;             
     P_image(ind_ght)=-p_not(j,ind_ght-idg+1);
     p_notg(j,:)=R0.*(p_not(j,:)./D1 + P_image(1:length(Time_Ax))./D2); % airgun signature
end