function [Signature]=agora_kernel(ArrayGeo)
 % Ozkan Sertlek, 2015, University of Leiden, 
 
%     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/>.

%     AGORA is a research algorithm which calculates the airgun array 
%     signatures to use in the sound mapping tools for the environmental 
%     risk assessment. The motivation of this algorithm is based on a PhD 
%     project (at Institute of Biology, University of Leiden) which aims to
%     model the soundscape of the Dutch North Sea. This project is funded 
%     by ZKO programme of NWO. The mathematical basement of the algorithm 
%     is described by "H.Ö.Sertlek and M.A. Ainslie, AGORA: Airgun source 
%     signature model: its application for the Dutch seismic surveys, 
%     Conference Proceedings of UAC 2015, Crete, Greece". 

%     The developers of AGORA accepts no responsibility or liability whatsoever
%     with regard to the results of AGORA if they are used in any commercial or scientific projects.
%     It is a research algorithm and still needs more advanced approach 
%     to increase its accuracy for entire frequency band. 
%
%     If you have any feedback or questions, please contact agora.model@gmail.com.  

% INPUT : ArrayGeo is a structure which provides input for the airgun geometry. 
% ArrayGeo.Pg= firing presssure [psi]
% ArrayGeo.Vg= airgun chamber size [in3]
% ArrayGeo.zg= firing depth [m]
% 
% ArrayGeo.N= Number of airguns
% ArrayGeo.x= x location of gun  [m]
% ArrayGeo.y= y location of gun [m]
% ArrayGeo.zr= receiver depth [m]
 
% OUTPUT: Signature is a structure which has the airgun signatures as 
% Signature.p_unperturbed % notional signature without bubble interactions
% Signature.Time_Ax    % time axis
% Signature.p_notional % notional signature with bubble interactions
% Signature.p_bi_ghost % notional signature with the surface ghost and bubble interactions

%% version notes
% AGORA version 2.0
% - The bugs related with the derivative of the bubble wall velocity are fixed in agora_kernel.m 
% - Signature calculation for International Airgun Modelling Workshop (IAMW) 2016 are added
% - The same version number from IAMW 2016 is used. 
%%
hwa=waitbar(0,'the simulation is started...')
input_to_AGORA                      %inputs for the simulations
%% read inputs from file 
env_in=load('environment.ain');

% calibration coefficients 
tao0=env_in(5);
beta=env_in(6);
eta=env_in(7);
k=env_in(8);

dt=env_in(9); % time resolution
cv=718*Ma;

% unit conversion
for gun_id=1:(ArrayGeo.N)  %loop over individual guns
%unit conversions to SI system
Pg=ArrayGeo.Pg(gun_id)*psi_to_Pa+P_atm;
Vg=ArrayGeo.Vg(gun_id)*in3_to_m3;
P_inf=P_atm +rhow_inf*g*ArrayGeo.zg(gun_id);

% intial parameters
ti=0;
zg0=ArrayGeo.zg(gun_id);% not disturbed gun depth
Tg0=Tw*(1+Pg/139e6);    % initial temperature in the gun chamber
%m0=P_inf*Vg/ (Rg*Tw);      
mg=Pg*Vg/ (Rg*Tg0);     % total mass in the gun chamber
a0=(3*Vg/(4*pi))^(1/3); % initial bubble radius- derived from the volume of gun chamber 
U0=0;                   %-(P_inf-Pg)/(rhow_inf*c_inf); % initial bubble wall velocity
mb0=P_inf*Vg/ (Rg*Tw);  % intiial mass of the bubble

% initial values in the stracture form
BubbleData.Ro(gun_id,1)=a0;
BubbleData.Uo(gun_id,1)=U0;
BubbleData.To(gun_id,1)=Tg0;
y=[mb0 Tw a0 U0];   % iteration start values
%% iterative solution of the differential equation system
    waitbar(((gun_id-1))/(ArrayGeo.N),hwa,[num2str(gun_id) '. airgun of ' num2str(ArrayGeo.N)  ' airguns is being calculated...'])

for t=0:dt:tmax
    hold on
    waitbar(((gun_id-1)*tmax/dt+(t)/dt)/(ArrayGeo.N.*length(0:dt:tmax)),hwa,[num2str(gun_id) '. airgun of ' num2str(ArrayGeo.N)  ' airguns is being calculated...'])
    ti=ti+1;        % itetation index
    mb=y(1);     Tb=y(2);    a =y(3);   u =y(4);     % start values
    BubbleData.Ro(gun_id,ti)=a;
    Vb=(4/3)*pi*(a^3);
    if ti==1
        Pa=P_inf;           % for the first time step pressure at bubble wall equals to P_inf  
    else
        Pa= mb*(Rg*Tb)/Vb;  % Pressure is found from PV=mRT formula 
    end
    BubbleData.P_bubblewall(gun_id,ti)=Pa;  % pressure at the bubble wall (matrix form)
    
    % calculate perturbed density and sound speed
    rhow=rhow_inf*(((Pa+B)/(P_inf+B))^(-1/n));
    c=c_inf*((Pa+B)/(P_inf+B))^((n-1)/(2*n));
    
    % calculate the differentials  
    dVb=4*pi*u*(a^2);  % variation in Volume
    dQb=k*4*pi*(a^2)*(Tb-Tw);    % variation in heat  
    mgg=abs(mg-mb);            
    Pgg=mgg*(Rg*Tg0)/Vg;
     if mb<=eta*mg  && Pgg>=Pa   % depends on airgun efficiency factor (eta). Some air remains in the chamber
         dmb=tao0*(Vg^beta)*sqrt((mgg/Vg)*(Pgg-Pa));            % variation in mass        
     else
         dmb=0;
     end

    dTb=(Rg*Tb*dmb-dQb-Pa*dVb)/(mb*cv);                          % variation in the temperature
    dH=( Rg/ (Vb^2))*(Tb*Vb*dmb + mb*Vb*dTb -mb*Tb*dVb)/(rhow);  % variation in the enthalpy
   
  
    H=n*(P_inf+B)*( ((Pa+B)/(P_inf+B))^((n-1)/n) -1)/(rhow_inf*(n-1));  % calculates the enthalpy
    Atemp=(1+u/c)*H + (a/c)*dH*(1-u/c) - 1.5*(u^2)*(1-u/(3*c)); 
    Btemp=a*(1-u/c);
    du(ti+1)=Atemp/Btemp;    % variation in the bubble wall velocity

    % second derivative
    d2u=(du(ti+1)-du(ti))/dt;
    
    DR=u*dt+du(ti+1)*(dt^2)/2;      % change in bubble radius
    DU=du(ti+1)*dt+d2u*(dt^2)/2;    % change in bubble velocity
    
    % calculate the new variables 
    BubbleData.Ro(gun_id,ti+1)=a+DR;
    BubbleData.Uo(gun_id,ti+1)=u+DU;
    BubbleData.To(gun_id,ti+1)=Tb+dTb*dt;
    BubbleData.mo(gun_id,ti+1)=mb+dmb*dt; 
    
    % variation in the source depth
    sumR=trapz(0:dt:(ti*dt),(BubbleData.Ro(gun_id,1:ti+1)).^3);
    dh=2*g*(BubbleData.Ro(gun_id,ti+1)^-3)*sumR*dt;
    zg=zg0+dh;
    % modified P_inf with the new source depth
    P_inf=P_atm +rhow*g*zg;
    
    % next iteration input
    y(4)=u+DU;
    y(3)=BubbleData.Ro(gun_id,ti+1);
    y(2)=Tb+dTb*dt;
    y(1)=mb+dmb*dt;
end
H=n*(P_inf+B).*(((BubbleData.P_bubblewall(gun_id,:)+B)./(P_inf+B)).^((n-1)/n) -1)./(rhow_inf*(n-1));

Signature.p_unperturbed(gun_id,:)=rhow_inf.*BubbleData.Ro(gun_id,1:end-1).*(H+(BubbleData.Uo(gun_id,1:end-1).^2)./2);
Signature.Time_Ax=1000*(0:dt:tmax); % ms
end
close(hwa)

%% bubble interactions 
[Signature.p_notional Signature.p_bi_ghost]=add_bubble_interactions_dist(BubbleData,Signature, ArrayGeo,cw,rhow_inf,dt,r);
     
%% FFT
%fmax=5000; % at the high frequencies the accuracy is not clear...
for j=1:ArrayGeo.N
    input_signature=Signature.p_notional(j,:);
    sig_len=length(input_signature); 
    t_pad=ceil(sig_len/10);  % applies zero padding
    input_signature((t_pad+1):sig_len+t_pad)=input_signature;
    input_signature(1:t_pad)=0;
    input_signature(2*sig_len)=0;
   
%     [Signature.F0i Signature.P_out(j,:)]=fft_func(input_signature,dt); % 
% for memory reduction in the directivity calculation
    fmax=3000; % at the high frequencies the accuracy is not clear...
    [f0 P_out_raw]=fft_func(input_signature,dt); % 
    F0id=find(f0<=fmax);
    P_out_raw=P_out_raw(F0id);
    F0=f0(F0id);
    Signature.F0i=1:F0(end);
    rP_out=interp1(F0,real(P_out_raw),Signature.F0i);
    iP_out=interp1(F0,imag(P_out_raw),Signature.F0i);
    Signature.P_out(j,:)=rP_out+1i*iP_out;
    clear P_out_raw rP_out iP_out
end
