calc_vel_int_Gz

 compute:
  \cG^z(\Delta,z) = 2\pi  \int k\id k\; J_0(k\Delta) \cG^z(k,z)

0001 function [I]=calc_vel_int_Gz(p);
0002 %
0003 %
0004 % compute:
0005 %  \cG^z(\Delta,z) = 2\pi  \int k\id k\; J_0(k\Delta) \cG^z(k,z)
0006 %
0007 %
0008 
0009 
0010 Nk=2000;
0011 k_vector=linspace(min(p.k),max(p.k),Nk);
0012 delta=p.delta_grid;
0013 t_out=1:800;
0014 I=zeros(length(p.w),length(t_out),length(delta));
0015  
0016 %the_i=39;
0017 for (i=1:length(p.w)); %%1:length(p.w))
0018   
0019   tic;
0020   the_w=p.w(i);
0021   
0022   %
0023   % compute the greens functions
0024   %
0025   %str=[[p.out_directory] 'greens/G2rs'],
0026   %a=load(str);
0027   %p=a.p;
0028   %G=a.G;
0029   %t_out=a.t_out;
0030   %k=p.k;
0031   %omega=p.w;
0032   omega=the_w;
0033   k=k_vector;
0034   SRC=1;       %% vertical momentum source
0035   OBS=1;       %% observe vertical velocity
0036   FIXED_OBS=1; %% fixed observation height
0037   F=0;         %% don't reload the eigenfunctions
0038   %% except for the first time
0039   if (i==1),F=1;end;
0040 
0041   t_s=p.t_src;
0042   t_obs=p.t_obs;
0043   USE_OTF=1;   %% this integral involves curly G, so use the OTF
0044   
0045   %
0046   % compute curly G
0047   %
0048   [t_out,G,c]=calc_green6(omega,k,USE_OTF,SRC,OBS,t_obs,FIXED_OBS,p,F);
0049   G=squeeze(G);
0050 
0051   
0052   %
0053   % multiply by k
0054   
0055   for (j=1:length(t_out))
0056       G(j,:)=G(j,:).*k_vector*1e-8;
0057   end
0058   
0059   if (find_min_width(G(100,:))<5)
0060     write_warning(p,'found a width of less than five pixels in calc_vel_int_Gz');
0061     write_warning(p,['at omega grid point of ' num2str(i)]);
0062   end;
0063   
0064  
0065   %
0066   % bessel transform
0067   %
0068   tmp=(1e-8)*uni_fast_bessel_trans2(0,k_vector,G,delta);
0069   I(i,:,:)=2*pi*reshape(tmp,1,length(t_out),length(delta));
0070   ttt=toc;
0071   if (noisy(p))
0072 disp(['calc_vel_int_Gz: i= ' num2str(i) ': time= ' num2str(ttt)]);
0073   end;
0074   
0075 end;
0076 
0077 
0078 %
0079 % save the output
0080 %
0081 
0082 
0083 str=p.out_directory;
0084 if (~exist(str))
0085   eval(['! mkdir ' str]);
0086 end;
0087 
0088 str=[p.out_directory 'integrals/'];
0089 if (~exist(str))
0090   eval(['! mkdir ' str]);
0091 end;
0092 
0093 w=p.w;
0094 save([p.out_directory 'integrals/vel_Gz'],'I','k_vector','w');
0095 save([p.out_directory 't_out'],'t_out','c');