test_horiz_int_v2

 interpolation method

0001 function [xg,offset_vector,sumh,sumv]=test_horiz_int_v2(p,the_dist,z_index);
0002 
0003 
0004 %
0005 % interpolation method
0006 %
0007 method='cubic';
0008 
0009 iter_calc_vel_int_Gz(p);
0010 iter_calc_vel_int_Iz(p);
0011 iter_calc_vel_int_Gh(p);
0012 iter_calc_vel_int_Ih(p);
0013 
0014 
0015 str=[p.out_directory 'integrals/vel_Iz'],
0016 a=load(str);
0017 Iz=a.I;
0018 dIz=a.dI;
0019 clear a;
0020 
0021 str=[p.out_directory 'integrals/vel_Gz'],
0022 a=load(str);
0023 Gz=a.I;
0024 t_out=a.t_out;
0025 clear a
0026 
0027 str=[p.out_directory 'integrals/vel_Ih'],
0028 a=load(str);
0029 Ih=a.I;
0030 dIh=a.dI;
0031 clear a;
0032 
0033 str=[p.out_directory 'integrals/vel_Gh'],
0034 a=load(str);
0035 Gh=a.I;
0036 clear a
0037 
0038 
0039 a=load([p.model_file]);
0040 m=a.m;
0041 z=interp1(m.t,m.r,t_out);
0042 clear a;
0043 
0044 %
0045 %
0046 % slice through the kernel for Vx
0047 %
0048 
0049 max_delta=max(p.delta);
0050 N=200;
0051 %M=10;
0052 xg=linspace(-N*max_delta,N*max_delta,N*p.N_gridx);
0053 length(xg),
0054 
0055 %
0056 %
0057 %
0058 Dist=p.delta(the_dist);
0059 
0060 
0061 %offset_vector=linspace(0,M*max(p.delta),M*p.Ngridy);
0062 offset_vector=linspace(0,0,1);
0063 length(offset_vector),
0064 
0065 sumh=zeros(size(offset_vector));
0066 sumv=zeros(size(offset_vector));
0067 
0068 for (wo=1:length(offset_vector))
0069   if (~rem(wo,100)),disp([num2str(round(100*wo/length(offset_vector))) '% done']);end;
0070   offset=offset_vector(wo);
0071   %
0072   % distance from one and two to r
0073   %
0074   %
0075 
0076   D1=sqrt(  (xg+Dist/2).^2 + offset^2 );
0077   D2=sqrt(  (xg-Dist/2).^2 + offset^2 );
0078 
0079   %
0080   % x and y components of \hat{\bdist}_1
0081   %
0082 
0083   D1x= (xg+Dist/2)./sqrt( (xg+Dist/2).^2 + offset^2);
0084   D1y= offset.* (sqrt( (xg+Dist/2).^2 + offset^2)).^(-1);
0085   
0086   D2x=(xg-Dist/2)./sqrt( (xg-Dist/2).^2 + offset^2);
0087   D2y=offset.* (sqrt( (xg-Dist/2).^2 + offset^2)).^(-1);
0088   
0089   D1xx=offset^2 ./ D1.^3;
0090   D2xx=offset^2 ./ D2.^3;
0091   D1yx=-offset*(xg+Dist/2)./D1.^3;
0092   D2yx=-offset*(xg-Dist/2)./D2.^3;
0093 
0094   D2yx(isnan(D2yx))=0;
0095   D2xx(isnan(D2xx))=0;
0096   D1xx(isnan(D1xx))=0;
0097   D1yx(isnan(D1yx))=0;
0098   D1x(isnan(D1x))=0;
0099   D2x(isnan(D2x))=0;
0100   D1y(isnan(D1y))=0;
0101   D2y(isnan(D2y))=0;
0102 
0103 
0104   d=p.delta_grid;
0105   w=p.w;
0106   if (length(w)>1)
0107     dw=w(2)-w(1);
0108   else
0109     dw=1;
0110   end;
0111   
0112   [Nw,Nt,Nd]=size(Iz);
0113 
0114 
0115   %
0116   % either loop over all frequencies, or just one
0117   %
0118   
0119   if (~isfield(p,'single_tw'))
0120     wlist=1:length(w);
0121   else
0122     wlist=p.single_tw;
0123   end;
0124 
0125 
0126   for (i=wlist)
0127 
0128     the_omega=w(i);
0129  
0130     Gz1=interp1(d,squeeze(Gz(i,z_index,:)),D1,method);Gz1(D1>=max(d))=0;
0131     Gz2=interp1(d,squeeze(Gz(i,z_index,:)),D2,method);Gz2(D2>=max(d))=0;
0132     Gh1=interp1(d,squeeze(Gh(i,z_index,:)),D1,method);Gh1(D1>=max(d))=0;
0133     Gh2=interp1(d,squeeze(Gh(i,z_index,:)),D2,method);Gh2(D2>=max(d))=0;
0134     
0135     
0136     Ih_at_omega=squeeze(Ih(i,z_index,:));
0137     Iz_at_omega=squeeze(Iz(i,z_index,:));
0138     dIz_at_omega=squeeze(dIz(i,z_index,:));
0139 %    dIz_at_omega=gradient(Iz_at_omega,d*1e8);
0140     
0141  %   dIh=gradient(Ih_at_omega,d);
0142    dIh_at_omega=squeeze(dIh(i,z_index,:));
0143     
0144     dIh1=interp1(d,dIh_at_omega,D1,method);dIh1(D1>=max(d))=0;
0145     Ih1=interp1(d,Ih_at_omega,D1,method);Ih1(D1>=max(d))=0;
0146     
0147     gi=dIh1.* D1x.^2 +Ih1.*D1xx;
0148     
0149     term1=sqrt(-1)*Gh2.*D2x.*gi;
0150     
0151    %
0152    %
0153    % second term
0154    %
0155    %   term2(j,:)=     sqrt(-1)*w(i)*Gh2(j,:).*D2y.*gradient(D1y.*Ih1(j,:),xg);
0156    %
0157      
0158 
0159   
0160 
0161 
0162    
0163    gi=dIh1.* D1x.*D1y+Ih1.*D1yx;
0164    
0165    term2=sqrt(-1)*Gh2.*D2y.*gi;
0166   
0167    
0168    %
0169    % third term
0170    %
0171    %
0172    %
0173    %term3(j,:)=    -sqrt(-1)*w(i)*Gz2(j,:).*gradient(Iz1(j,:),xg);
0174    
0175  %  dg=gradient(Iz_at_omega,d);
0176    deriv_g=interp1(d,dIz_at_omega,D1,method).*D1x;deriv_g(D1>=max(d))=0;
0177    
0178     
0179    term3=-sqrt(-1)*Gz2.*deriv_g;
0180    
0181    
0182    
0183    
0184    %
0185    %
0186    % first term cc
0187    %
0188    %
0189     
0190    dIh2=interp1(d,dIh_at_omega,D2,method);dIh2(D2>=max(d))=0;
0191    Ih2=interp1(d,Ih_at_omega,D2,method);Ih2(D2>=max(d))=0;
0192    gi=dIh2.* D2x.^2+Ih2.*D2xx;
0193    term1cc=conj(sqrt(-1)*Gh1.*D1x.*gi);
0194   
0195    %
0196    % second term cc
0197    %
0198    %
0199   
0200    gi=dIh2.* D2x.*D2y+Ih2.*D2yx;
0201   
0202    term2cc=conj(sqrt(-1)*Gh1.*D1y.*gi);
0203   
0204 
0205   
0206    %
0207    % third term cc
0208    %
0209    %    term3cc(j,:)=conj( -sqrt(-1)*w(i)*Gz1(j,:).*gradient(Iz2(j,:),xg) );
0210    %
0211   
0212    deriv_g=interp1(d,dIz_at_omega,D2,method).*D2x;deriv_g(D2>=max(d))=0;
0213    
0214    term3cc=conj(-sqrt(-1)*Gz1.*deriv_g);
0215 
0216 %   disp('**** std of the three terms:');
0217 %   term1(isnan(term1))=0;
0218 %   term2(isnan(term2))=0;
0219 %   term3(isnan(term3))=0;
0220   
0221 %   std(abs(term1)),
0222 %   std(abs(term2)),
0223 %   std(abs(term3)),
0224    
0225  %  term1cc(isnan(term1cc))=0;
0226  %  term2cc(isnan(term2cc))=0;
0227  %  term3cc(isnan(term3cc))=0;
0228    
0229  %  std(abs(term1cc)),
0230  %  std(abs(term2cc)),
0231  %  std(abs(term3cc)),
0232    
0233    
0234    
0235    vh=-2*1e-8*the_omega*(term1+term1cc+term2+term2cc);
0236    vh(isnan(vh))=0;
0237    
0238   
0239    vz=-2*the_omega*(term3+term3cc);
0240    vz(isnan(vz))=0;
0241 
0242    
0243   
0244    
0245    %disp(['std vh and vz are ' num2str(std(abs(vh(:)))) ' and ' num2str(std(abs(vz(:)))) ]);
0246    %disp(['sum(vh) and sum(vz) are' num2str(sum(vh)) ' and '  num2str(sum(vz)) ]);
0247    
0248  
0249  sumh(wo)=sum(vh);
0250  sumv(wo)=sum(vz);
0251    
0252   end;
0253 
0254 end;
0255 
0256 c=interp1(m.t,m.c,t_out);
0257 rho=interp1(m.t,m.rho,t_out);
0258 
0259   
0260 %dCc=-2e-8*(term3+term3cc);
0261 
0262 dx=xg(2)-xg(1);
0263 if (length(offset_vector)>1)
0264   dy=offset_vector(2)-offset_vector(1);
0265 else
0266   dy=1;
0267 end;
0268 
0269 scale=rho(z_index)*(2*pi)^7*calc_mzero(p.w,p);
0270 sumh=scale * sumh;
0271 sumv=scale * sumv;
0272 
0273 %dCc(isnan(dCc))=0;
0274 close all;plot(xg,imag(vz),'+');
0275 %disp(['total integral is ' num2str(2*sum(sumh+sumv)*dx*dy*1e16)]);
0276 disp([' h integral is ' num2str( (2*sum(sumh)+sumh(1))*dx*dy*1e16)]);
0277 disp([' z integral is ' num2str( (2*sum(sumv)+sumv(1))*dx*dy*1e16)]);
0278