function [F_delta,pos,freqp,df,times]=seismic_response_primary_list(pathNC, ...
    dpt, lat,lon, CgR,  Q, date1, date2, dt, lono, lato, var, deep)

%pathNC=pathref1;
%dpt=dpt1; 
%var=varef;
%deep=4000.;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This functions computes synthetic seismic spectra due to primary mechanism
% WITH SOME ASSUMPTIONS ON THE DIRECTIONAL SPECTRUM... 
% at the point of longitude lono and latitude lato
%
% Last modified: 29 Sep 2013
%
% History     : created September 2013 by F. Ardhuin
%
% Input parameters: 
% - pathNC    : file path where the *.p2f NetCDF files are stored
% - dpt       : array of water depths, should have the same size as the
%               wave data
% - Q         : quality factor (function of frequency)
% - date1     : start date in matlab form =datenum(iy,im,id,ih,imin,isec)
% - date2     : end date in matlab form   =datenum(iy,im,id,ih,imin,isec)
% - dt        : time step of computation, in fraction of day
% - lono      : longitude of receiving station, where spectrum is estimated
% - lato      : latitude of receiving station, where spectrum is estimated
% - nmodes    : number of modes in calculation 
% - nmodesout : number of modes in output
%
% Output parameters: 
% - F_delta   : Spectral density of surface ground displacement in m^2/Hz
%               normalized so that sum(Ef*df) is the ground displacement
%               variance
% - pos       : centroid of sources 
% - freq      : seismic frequencies
% - df        : frequency increment
% - times     : matlab dates at which the spectrum is computed
%
% Method:
% Uses Hasselmann (1963): noise source for a straight coastline with 
% constant bottom slope (defined by Hp). 
% This function reads fr
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 1. Defines constants and coefficients
%
g    =9.81;     % gravity          in m/s^2
rhow =1026;     % density of water in kg/m^3
rhos =2600;     % density of rock  in kg/m^3
betas=2800;     % shear wave velocity of rock in m/s
R_E=4E7/(2*pi); % Earth radius     in m
pos=0.;
lg10=log(10);
%
%
% NB: the F_p2(K,f) is computed from E(f) in the NetCDF files 
%     spectral density of surface pressure (in meters)
%     such that (rhow*g)^2*F_p2(K,f)*df*dKx*dKy is
%     the variance of the pressure (in Pa^2) with df the 
%     ocean wave frequency interval.
%
%
% factor1 is the constant part of the integrand, 
% factor1=W^2*2*pi^2*(g*rhow)^2/(c^2*sigma)*sin(alpha)
%
factor1=pi*(g*rhow/rhos)^2/(betas^5*R_E); 
%
%
% This source is attenuated during propagation 
% 


%
% 2. Reads netcdf file to get dimensions ... .
%
nt=floor((date2-date1)/dt)+1;       % number of time steps
times=zeros(nt,1);
date=datevec(date1);
filename=sprintf('ww3.%.4d%.2d_ef.nc',date(1),mod(date(2),100))
fullfile=[pathNC filename]
[lat,lon,time,freqp,P2f,unit,scale,fillv,MAPSTA] ...
         =read_WWNCf_noscale(fullfile,var,date1,1);
     
%
nf=length(freqp);
Hp=0.01;  % slope ... 
omega=(2*pi).*freqp;                     % seismic radian frequency omega=2*2*pi/T


%
% Tabulation of Ef with scale factor
%
unit0=unit;
scale0=scale;
fillv0=fillv;
imin=-2^15+1;
imax=2^15; 
tabEf=linspace(imin,imax,imax-imin+1);
switch unit 
   case 'log10(m2s+1E-12)'
      tabEf=exp(lg10.*tabEf.*scale)-(1E-12-1E-16);
   case  'log10(m2 s+1E-12)'
      tabEf=exp(lg10.*tabEf.*scale)-(1E-12-1E-16);
   otherwise 
       warning(['Unexpected units for elevation spectrum: ' unit]);
end
%
tabEf(fillv+imax)=NaN;
%
%
%
nx=size(dpt,1);
ny=size(dpt,2);
la2=repmat(lat',nx,1);
lo2=repmat(lon,1,ny);
%
% distances from station
%
alpha=dist_sphere(lono,lo2,lato,la2).*(pi/180);  % spherical distance
coslat=repmat(cos(lat*pi/180)',size(lon),1);     % cos(latitude): 2D array
coslat1=reshape(repmat(cos(lat*pi/180)',size(lon),1),nx*ny,1);     % cos(latitude): 1D
sinlat1=reshape(repmat(sin(lat*pi/180)',size(lon),1),nx*ny,1);     % sin(latitude): 1D
coslon1=reshape(repmat(cos(lon*pi/180),1,ny),nx*ny,1);
sinlon1=reshape(repmat(sin(lon*pi/180),1,ny),nx*ny,1);

la2=repmat(lat',nx,1);
lo2=repmat(lon,1,ny);

%
% Elementary surface
%
dlon=(lon(nx)-lon(1))/(nx-1);
dlat=(lat(ny)-lat(1))/(ny-1);
coslat=repmat(cos(lat*pi/180)',size(lon),1);     % cos(latitude): 2D array
dA=R_E.^2.*coslat.*(dlon*dlat*(pi/180).^2);    % Area of surface element
nd=36;
theta=90-linspace(0,360*(nd-1)./nd,nd);

%
% 3. Finds the sources of primary noise as the land-sea boundary 
%
% 3.a resolved land  
  [bound,dirs,type]=seismic_primary_find_sources(dpt,nd);
  indb1=find(bound==1);
  nbound1=length(indb1);
%
% 3.b unresolved islands
%
  [Sx1,Sy1]  = read_obstr('glob_30m.obst',nx,ny);
  obst=(Sx1+Sy1)'/200;
  indo=find(obst>0 & bound ==0 & isfinite(dpt)==1);
  nobst=length(indo);
  nbound=nbound1+nd*nobst;
  ipd=zeros(nbound,1);
  indb=zeros(nbound,1);
  dkcross=zeros(nbound,1);
  dptp=zeros(nbound,1);

  dtor=pi/180;

%
% Fills array index for resolved land
%
  dlon=lon(2)-lon(1);
  dlat=lat(2)-lat(1);
  for ind=1:nbound1;
      ipd(ind)=1+mod(round(dirs(indb1(ind)))-1,nd);
      sind=abs(sin(theta(ipd(ind))*dtor));
      cosd=abs(cos(theta(ipd(ind))*dtor));
      dkcross(ind)=2*pi./((R_E*pi/180)*(dlon.*coslat(indb1(ind))*sind+dlat.*cosd));
      indb(ind)=indb1(ind);
      dptp(ind)=dpt(indb1(ind));

      % this is the wavenumber increment in cross-shore direction
      % corresponding to 2*pi / dl with dl the cell size. 

  end
  
%
% Fills array index for unresolved islands
%
  drad=pi/360;
  for i=1:nobst;
      for ith=1:nd
         ind=nbound1+(i-1)*nd+ith;
         indb(ind)=indo(i);
         % shifts point in direction theta(ith) ... 
         % looks for point output closest to shifted point
         ipd(ind)=ith;
         dptp(ind)=dpt(indo(i));
         
         sind=abs(sin(theta(ith)*dtor));
         cosd=abs(cos(theta(ith)*dtor));
         
         dkcross(ind)=2*pi./((R_E*pi/180)*(dlon.*coslat(indb(ind))*sind+dlat.*cosd)).*(obst(indb(ind))*pi/nd);  
      % the (obst(indb(ind)*pi/nd) factor corresponds to the reduction in shoreline
      % length compared to one full square... this assumes a round subgrid
      % island.
      end
  end
  
nf=length(freqp);
Hp=0.01;  % slope ... 
omega=(2*pi).*freqp;                     % seismic radian frequency omega=2*2*pi/T


%
% distances from station
%
alpha=dist_sphere(lono,lo2,lato,la2).*(pi/180);  % spherical distance

EftoP2f=zeros(nbound,nf);

% tabulation of dispersion relation 
Y1=0.02;
Y2=3.55;
dY=0.001;

Ytab=linspace(Y1,Y2,round((Y2-Y1)./dY)+1);
frY=sqrt(Ytab.*(g/100))./(2.*pi);
Xtab=dispNewtonTH(frY,100.).*100.;


for j=1:nf
   wnumdeep=omega(j).^2/g;
   kd=dispNewtonTH(freqp(j),deep);
   Cgd=omega(j)/kd*(0.5+kd*deep/sinh(2*kd*deep));
      Y=omega(j).^2.*dptp./g;
      X=Y;
      I1=find(Y <= 0.02);
      I2=find(Y > 0.02 & Y < 3.55);
      I3=find(Y > 3.55);
      % X is such that Y = X*tanh(X) 
      X(I1)=sqrt(Y(I1));
      X(I2)=Xtab(round((Y(I2)-Y1)./dY)+1);
      kl=X./dptp;
      Cgl=omega(j)./kl.*(0.5+kd*dptp./sinh(2.*kl.*dptp));
      % This corrects for the shoaling effect of waves between the 'deep' depth and 
      % the local depth. See Ardhuin et al. (Ocean Modelling 2014)
      facshoal=kd.*Cgl./(kl.*Cgd);
      if j==11
          facshoalj=facshoal;
      end
      EftoP2f(:,j)=2.*((g.^2./omega(j)^4)*(Hp/(64*pi)).*dkcross/kd)./(2*pi).*facshoal;
end



%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Defines seismic response function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Version 1: fit to the curves shown in Kedar et al. 2008
%            based on the tabulated values in Longuet-Higgins 1950
%
nC2=601;
x0=linspace(0,6,nC2);   % discretized values of sighoverbeta*100
                        % i.e. sig*h/beta goes from zero to 6
a=1.;
b=10;
c=0.4;
c1=0.6*(a+c.*x0)./(a+b.*(x0-0.85).^2)+0.28./(2.+1.*x0)-0.008*(x0-2.5);

a=6;
b=30;
c=0.5;
d=10;
c2=(0.30*(a+c.*x0)./(a+b.*(x0-2.77).^2)+0.40./(1.+1.*x0)+0.006*(x0-5)).*max(tanh((x0-1)*d),0) ...
    + 0.02.*exp(-(1.5*(x0-3.6)).^2)-max(0.02*(x0-5),0);

C2=c1.^2+c2.^2;     % first 2 modes: this is c^2
c0=sqrt(c1.^2);     % first mode only


readc=0;
if readc==1  
% Version 2
% Reads in the response functions c=sqrt(c1.^2+c2.^2+ ...)
% computed with sismo_stoneley_all.m
%
   nmodes = 1;
   all_c=load('Rayleigh_source.txt');
   taille=size(all_c);
   nC=taille(1);
   nmodein=taille(2)-1;
   if nmodes > nmodein
       nmodes=nmodein;
   end
   sighoverbetain=pi*all_c(:,1);
   nC2=1001;
   x0=linspace(0,10,nC2);   % discretized values of sighoverbeta*100
                            % i.e. sig*h/beta goes from zero to 6
   C2=x0.*0;
   for i=1:nmodes
       c=interp1(sighoverbetain,all_c(:,i+1),x0);
       I=find(isfinite(c));
       C2(I)=C2(I)+c(I).^2;
   end
   C2(1)=C2(2);
end

%
% Defines group speeds
%
if CgR < 0 
   all_Cg=load('Rayleigh_Cg.txt');
   taille=size(all_Cg);
   nCg=taille(1);
   nmodein=taille(2)-1;
   sighoverbetain=pi*all_Cg(:,1);
   nC2=1001;
   x0=linspace(0,10,nC2);   % discretized values of sighoverbeta*100
                            % i.e. sig*h/beta goes from zero to 6

   CgRcmax=x0.*0;
   CgRcmax=interp1(sighoverbetain,all_Cg(:,2),x0);
   c0=interp1(sighoverbetain,all_c(:,2),x0);
   for i=2:nmodes
       c=interp1(sighoverbetain,all_c(:,i+1),x0);
       Cgi=interp1(sighoverbetain,all_Cg(:,i+1),x0);
       I=find(isfinite(c) & isfinite(c0));
       I2=find(c(I) > c0(I) );  
       CgRcmax(I(I2))=Cgi(I(I2));
       c0(I(I2))=c(I(I2));
   end
   I=find(isfinite(CgRcmax) == 0);
   CgRcmax(I)=CgRcmax(I(1)-1);
else
    CgRcmax=zeros(nC2)+CgR;
end

%figure(1)
%plot(x0,CgRcmax,'k',sighoverbetain,all_Cg(:,2),'r--')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Uses wave model output and computes the seismic response
%

ntrmax=8;  % Number of time steps read at once (reduced I/O time)
itr=0;
for it=1:nt
    now=date1+(it-1)*dt;
    date=datevec(now);

%
% Reads NetCDF ef files 
%
    filename=sprintf('ww3.%.4d%.2d_ef.nc',date(1),mod(date(2),100));
%
% Reads maps of sources (frequency-integrated)
%
    fullname=[pathNC filename];
    yname=sprintf('%.4d',date(1));
    iy=date(1);
    if (mod(it,ntrmax)==1)
      ntr=min(ntrmax,nt-it+1);  
      [lat,lon,time,freqp,Eft,unit,scale,fillv,MAPSTA] ...
         =read_WWNCf_noscale(fullname,var,now,ntr);

      msg=sprintf(': %d out of %d ', it, nt);
      fprintf('reading NetCDF time step %s, in file %s \n',msg,fullname)
      itr=1;
      ntr=size(Eft,4);
    else
      itr=itr+1;
      %[it itr]
    end
    times(it)=now;
    if unit == unit0 & scale == scale0 & fillv == fillv0 
       Ef=tabEf(imax+int32(squeeze(Eft(:,:,:,min(itr,ntr)))));
    else
       'ERROR: CHANGE IN SCALE FACTOR OR UNITS IN NEW *.nc FILE'
    end
%

% pcolor(lon,lat,double(squeeze(Ef(:,:,10)))');shading flat;colorbar;

% 
%
% Computes amplification and attenuation for the first time step
%
    if it==1
        nf=length(freqp);
        pos_noref=zeros(nf,nt,4);
        coeff=zeros(nbound,nf);
        attenuation=zeros(nbound,nf);
        
        for j=1:nf
            omehoverbeta=0.*omega(j).*dpt(indb)./betas;   % set to zero because interaction in shallow water
            i=1+int16(100.*omehoverbeta);                % gets the nearest discretized omega*h/beta
            I=find(i>nC2);
            i(I)=nC2;
            coeff(:,j)=(factor1.*omega(j).*C2(i))./sin(alpha(indb)');
            attenuation(:,j)=exp(-1.*omega(j).*alpha(indb).*(R_E./(abs(CgR)*Q(j))));

            %imode = 1;
            %c=interp1(sighoverbetain,all_c(:,imode+1),x0).^2./C2;
            %CgRm=x0.*0+abs(CgR);
            
               
               I2=find(isfinite(c)==0 | C2 == 0);
               c(I2)=0;
        end
        alpha1D=reshape(alpha(indb),nbound,1);   % angular distance between source and station
        dA3D=repmat(dA(indb),[1,nf]);         % elementary Earth area
        coeff_all=coeff.*attenuation.*dA3D;
    end 
    P2f = zeros(nbound,nf);
    mat1D=zeros(nbound,1);
    for j=1:nf
       Ef2=Ef(:,:,j);
       P2f(:,j) = Ef2(indb).*EftoP2f(:,j);
    end

    source=coeff_all(:,:).*P2f(:,:);             
%
% Only takes points that are in the ocean (source is finite) and not too
% close (singularity for alpha=0) 
%
    I=find(isfinite(squeeze(source(:,1))) & sin(alpha1D) > 0.002);
    F_delta(:,it)= sum(source(I,:),1);
    if it==1
       figure(102)
       j=11;
       set(gca,'FontSize',16);
       title('P2f from ef');
        scatter(lo2(indb),la2(indb),20,P2f(:,j),'LineWidth',2);
        save 'ef_P2f' indb lo2 la2 P2f source coeff_all EftoP2f dkcross facshoalj
    end
end
nf=length(freqp);
xfr=exp(log(freqp(nf)/freqp(1))/(nf-1));  % determines the xfr geometric progression factor
df=freqp*0.5*(xfr-1/xfr);% frequency intervals in wave model times 2


%end

