calc_green9

0001 function [g]=calc_green9(omega_out,k_out,USE_OTF,SRC,OBS,t_obs,t_src,p,FORCE_RELOAD)
0002 %
0003 %
0004 %
0005 %
0006 %
0007 %
0008 %
0009 %
0010 %
0011 
0012 
0013 
0014 persistent already_loaded V W omega k model t c
0015 
0016 
0017 if (~length(already_loaded) | FORCE_RELOAD )
0018   %
0019   % this only runs the first time you call this routine
0020   % or if you explicitly ask for it
0021   % load the eigenfunctions and normalizes them
0022   %
0023   if noisy(p),disp('loading');end;
0024 
0025 
0026 
0027   already_loaded=1;
0028   %
0029   % this loads the eigenfunctions
0030   %
0031   % the vector velocity eigenfunction is
0032   %  \vel = W\hat{z}+
0033   [V,W,omega,k,t,model]=load_eigs2(p);
0034 
0035   if noisy(p),disp('done loading');end;
0036   
0037   
0038   c=interp1(model.t,model.c,t);
0039 end;
0040 
0041 
0042 %
0043 % have eigefunctions for this k
0044 %
0045 I=find( k_out > min(k) & k_out < max(k));
0046 
0047 
0048 %
0049 % get the size of eigenfunction set
0050 %
0051 [MAX_N,Nk,Nt]=size(V);
0052 
0053 %
0054 % get the source and receiver vectors all set up
0055 %
0056 
0057 
0058 if (SRC==0)
0059   S=-sqrt(-1)*my_interp(t,V,t_src,0,1,3);
0060 elseif(SRC==1)
0061   S=my_interp(t,W,t_src,0,1,3);
0062 elseif(SRC==2)
0063   S=my_interp(t,divergence2(k,V,W,t,c),t_src,0,1,3);
0064 elseif(SRC==3)
0065   S=my_interp(t,W,t_src,1,-c.^(-1),3);
0066 else
0067   error(['dont know how to do SRC=' num2str(SRC)]);
0068 end;
0069 
0070 
0071 if (OBS==0)
0072   Ob=my_interp_v2(t,V,t_obs,3);
0073 elseif(OBS==1)
0074   Ob=my_interp_v2(t,W,t_obs,3);
0075 elseif(OBS==2)
0076   Ob2=divergence2(k,V,W,t,c);
0077   Ob=my_interp_v2(t,Ob2,t_obs,3);
0078 elseif (OBS==3)
0079   Ob=my_interp(t,W,t_obs,1,-c.^(-1),3);
0080 elseif (OBS==4)
0081   Ob=my_interp(t,V,t_obs,1,-c.^(-1),3);
0082  else
0083   error(['dont know how to do OBS=' num2str(OBS)]);
0084 end
0085 
0086 %
0087 %
0088 % now interpolate onto the grid in k that we want the answer on
0089 %
0090 
0091 tmp_S=zeros(MAX_N,length(k_out));
0092 tmp_S(:,I)=fast_interp2(k,S,k_out(I));
0093 S=tmp_S;
0094 
0095 tmp_Ob=zeros(MAX_N,length(k_out));
0096 tmp_Ob(:,I)=fast_interp2(k,Ob,k_out(I));
0097 Ob=tmp_Ob;
0098 
0099 
0100 omega_int=zeros(MAX_N,length(k_out));
0101 omega_int(:,I)=my_interp(k,omega,k_out(I),0,1,2);
0102 if (~isfield(p,'quiet')),
0103   disp('done with interpolation');
0104 end;
0105 
0106 %
0107 % allocate space for the result
0108 %
0109 g=zeros(length(omega_out),length(k_out));
0110 
0111 
0112 
0113 
0114 %close all;
0115 %subplot(1,3,1);imagesc(abs(Ob));
0116 %subplot(1,3,2);imagesc(abs(S));
0117 %subplot(1,3,3);imagesc(abs(omega_int));
0118 
0119 
0120 for j=1:MAX_N
0121   for i=I
0122     A=Ob(j,i)*S(j,i);
0123     w=omega_int(j,i);
0124     A=reshape(A,1,length(A),1);
0125     if real(w)>p.max_mode_frequency*2*pi/1e3
0126       A=0;
0127     end
0128     
0129     for (l=1:length(omega_out))
0130       
0131       tmp=1/(2)*(  -conj(A)*(omega_out(l)+conj(w))^(-1)/conj(w) + A*(omega_out(l)-w).^(-1)/w );
0132       g(l,i)=g(l,i)-tmp; %reshape(t,1,length(t),1);
0133       
0134     end
0135   end
0136 end
0137 
0138 
0139 %
0140 % the OTF correction
0141 %
0142 if USE_OTF
0143   
0144   
0145   OTF=get_OTF(k_out,p.which_otf);
0146  
0147   for l=1:length(omega_out)
0148       g(l,:)=g(l,:).*reshape(OTF,1,length(OTF));
0149   end;
0150   
0151   %%
0152   %% then the data analysis filter
0153   %%
0154 
0155   F=make_filter(omega_out,k_out,p).';
0156    
0157   g=g.*F;
0158 
0159 end;
0160 
0161 
0162 
0163 g=g/((2*pi)^3);