k-Wave User Forum » Topic: Simulation "layer by layer"

k-Wave User Forum » Topic: Simulation "layer by layer" http://www.k-wave.org/forum/topic/simulation-layer-by-layer Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:39:03 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer/page/2#post-7576 Sat, 06 Jun 2020 21:44:52 +0000 Bradley Treeby 7576@http://www.k-wave.org/forum/ In case it's helpful, I dug out some of my old test codes from when Anthony was working on this and tidied one up. You can download it <a href="http://www.k-wave.org/downloads/dirichlet_divided_simulation_3D_cw.m">here</a>. This particular example is for a linear CW simulation with a focused transducer, and gives errors on the order of a few percent for a single grid point source plane. There is also a slightly different example in the <a href="http://www.k-wave.org/documentation/example_tvsp_equivalent_source_holography.php">Equivalent Source Holography Example</a>. DLam on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer/page/2#post-7540 Tue, 02 Jun 2020 19:34:16 +0000 DLam 7540@http://www.k-wave.org/forum/ Thanks for your experience and advice, Anthony. Yes, from my own work I noticed that changing the source plane thickness doesn't seem to change the peak intensity in a predictable way. I will try CFL = 0.5 next, and using the acoustic holography method (calculateMassSource set of functions?) to generate the source planes. Anthony on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer/page/2#post-7538 Tue, 02 Jun 2020 08:57:02 +0000 Anthony 7538@http://www.k-wave.org/forum/ Hi, Unfortunately, I did not find a general solution on my side, and did not have the time to dig further. From what I remember, I tweaked the parameters until it gave me appropriate results for my source. In particular, if I remember well, the errors were not decreasing with decreasing CFL (in an example cite above in the discussion, results were better with CFL = 0.5 than 0.3), and it was not monotonous either with increasing thickness of source "plane". I don't remember, but there might be something related to the staggered grid though, (see user manual for details and maybe on the forum, I think we discussed something on this matter somewhere)... It needs to be accounted for properly. Otherwise, you may try the equivalent source described by Brad in the paper cited above... Sorry not to be of more help :-( Maybe Brad's team can help. DLam on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer/page/2#post-7533 Thu, 28 May 2020 19:52:20 +0000 DLam 7533@http://www.k-wave.org/forum/ Hi Brad, Anthony, I've been doing multi-stage simulations as well. I have compared my multi-stage simulation (MULTI) results to the results obtained from a 'control' single-stage simulation (SINGLE), but I’m observing something odd in my results. For me, I am struggling with the transition between stages. The result is that the peak pressure is about 12% higher in MULTI than in SINGLE, as well as being offset in space (for MULTI, the focus is closer to the transducer surface than for SINGLE). The simulation was performed in a nonlinear homogeneous acoustic medium. The single-stage simulation was performed at a spatial resolution of 8.2 PPW, with CFL around 0.3. For the multi-stage simulation, all stages had the spatial resolution of 8.2 points per wavelength, but each stage was decreased in size as the acoustic beam narrowed. CFL was also around 0.3. Similar to Anthony, I recorded the source planes from the last several rows of a stage, then used them as the source in the following stage (the white bars). Source planes were 9 gridpoints thick, and were placed 10 gridpoints away from the non-PML domain boundary (so if the PML is 20 gridpoints thick, the source plane is 30 gridpoints away from the computational boundary edge). I record the last 4 cycles, and use interpftn to interpolate the recorded sensor pressures into becoming the new source plane. I looked at reducing the CFL number from 0.3 to 0.1 in MULTI. Peak intensity is closer to results obtained from SINGLE, but still too high (by approx 10% and the MULTI focus is still closer to the transducer surface than SINGLE). I am wondering if the interpolation from sensor to source is wrong somehow? Should I be using a windowing function, or time-shifting the interpolation somehow? Daniel Anthony on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer/page/2#post-7474 Mon, 04 May 2020 08:08:08 +0000 Anthony 7474@http://www.k-wave.org/forum/ Hi Brad, thank you for your answer, it is getting clearer in my mind. Warm regards, Anthony Bradley Treeby on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer/page/2#post-7459 Fri, 01 May 2020 08:51:18 +0000 Bradley Treeby 7459@http://www.k-wave.org/forum/ Hi Anthony, The errors are because the Dirichlet replacement in k-Wave is not exactly setting a boundary condition in a formal sense, just replacing the pressure values at discrete time points. If you run a simulation in 1D and record the pressure everywhere, then re-emit the recorded field from a single point using a Dirichlet source and record the pressure everywhere, and the compare the spatial gradients of the two fields at different time snapshots, you'll see they're not quite right. We have an idea how to fix this by formally deriving a boundary value propagator (similar to how we have derived solutions for initial value problems and single frequency problems). Practically, your work around of applying the source term over a plane several grid points thick helps to alleviate this. The key advantage of an additive source is that it doesn't act like a boundary, so if you have a heterogeneous medium, the scattered waves are not reflected from the source plane (they are if you use a Dirichlet source). For the equivalent source, I've not really played around with the minimum distance between the "measurement" plane and the source plane. In principle it should work even if they're separated by a single grid point, although I'd advise trying it out as there might be strange behaviour if you're this close. Brad. jackYANG on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer/page/2#post-7433 Wed, 22 Apr 2020 09:35:23 +0000 jackYANG 7433@http://www.k-wave.org/forum/ Hi Anthony, Thank you for your reply. I have downloaded the article you mentioned. I will try my best to solve this problem,If there is any progress, I will contact you again. Anthony on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer/page/2#post-7432 Wed, 22 Apr 2020 08:41:21 +0000 Anthony 7432@http://www.k-wave.org/forum/ Hi jackYANG, Unfortunately, I did not find a general solution... But there is a simple (yet ugly) trick: you can make your "layers" overlap enough so the oscillations at the beginning of the layer L+1 vanishes before the z value corresponding to the end of the layer L :-) Then, when you merge all the layers in the zone where L and L+1 overlap, take the pressure values coming from the layer L where it overlaps witj L+1... See figure 2 of <a href="https://jtultrasound.biomedcentral.com/articles/10.1186/s40349-016-0056-9" rel="nofollow">https://jtultrasound.biomedcentral.com/articles/10.1186/s40349-016-0056-9</a>. More elegantly, maybe that using the equivalent source computation method described in Brad's paper (see his last post in this very conversation) solves the problem. Please tell my if you try it and it works ;-) I would be very interested! Hope it helps, Anthony jackYANG on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer/page/2#post-7431 Wed, 22 Apr 2020 01:56:05 +0000 jackYANG 7431@http://www.k-wave.org/forum/ Hi Anthony, I'm facing the same problem you had before ，just like ，I let z denote my main propagation axis. When I compute the integral of Iz (intensity along z) over each xy-plane , and plot this quantity versus z, I get very strong oscillations close to the source plane.How did you solve this problem?Can I add CFL to solve it? Anthony on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer/page/2#post-7427 Tue, 21 Apr 2020 08:32:28 +0000 Anthony 7427@http://www.k-wave.org/forum/ Hi Brad, I just read the article you mentioned and I must admit that I am not sure to understand why directly replacing the pressure values in the domain "leads to errors in the imposed spatial gradients" (although I don't doubt it, since the example you give is very convincing). I am especially curious about the example of the simulated field, where measurement errors cannot be responsible for these errors on gradients, since you set the "true" pressure values... By the way, are you talking about inconsistencies in the 3D spatial gradients (including voxels out of the measurement plane) or only inside this 2D plane? Then, about the interpretation of the method, I understand that it is "restoring" the spatial consistency by allowing some (little) changes in the plane where the measurements were performed to make them physically consistent with the spatial gradients as they should be according to the equations, but how far do you need to put your equivalent source from the measurement plane? Is there a minimal distance? Finally, in the example of the section III B, how does the "additive" mode in k-Wave behave? Does it give the same errors as the "Dirichlet" mode, since it seems a more "permissive" way to set the source term. Warmest regards to everyone in this troubled times, Anthony Bradley Treeby on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer/page/2#post-7097 Fri, 18 Oct 2019 20:33:15 +0000 Bradley Treeby 7097@http://www.k-wave.org/forum/ Hi umitarabul, We've written something on it <a href="http://bug.medphys.ucl.ac.uk/papers/2018-Treeby-IEEETUFFC.pdf">here</a>, but haven't tried it for elastic waves. Do you want to split your simulation within the solid? Brad. umitarabul on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer/page/2#post-7080 Fri, 18 Oct 2019 15:48:27 +0000 umitarabul 7080@http://www.k-wave.org/forum/ I am interested in such a multi-step approach in 2-D for elastic waves. Is there an update about the amplitude mismatch issue mentioned in this topic? Anthony on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer#post-5202 Tue, 18 Aug 2015 07:23:10 +0000 Anthony 5202@http://www.k-wave.org/forum/ Thanks for the tips! Indeed, it's afordable :-) Good luck too! Anthony Eyal on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer#post-5201 Mon, 17 Aug 2015 18:55:58 +0000 Eyal 5201@http://www.k-wave.org/forum/ Hi Anthony, Thanks for the info, I decided to run the memory intensive steps on an amazon cloud server, they are not that expensive these days - you can get a server with 244GB of memory for about 3.5$ an hour. Good luck :) Eyal Anthony on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer#post-5197 Sun, 09 Aug 2015 15:04:59 +0000 Anthony 5197@http://www.k-wave.org/forum/ Hi Eyal, Unfortunately, it increases a lot the memory requirements during the preprocessing phase as the size of the source term is Ns*Nt, with Ns the number of source points and Nt the number of time points. In 3D, with fine grid steps, it was impractical on my hardware to use this method. It would become possible if the source term could be set as one period and automatically repeated by k-wave until the end of the temporal grid (with a 'continuous' flag)... (*** This is subliminal feature request ***). However, I don't see why it would'nt work in the nonlinear case... Finally, I increased the CFL to 0.5 and it appeared to be better (not perfect but good enough for my purpose). The development team will probably have more explanations about all this after the holydays :-) and I will post it here if I find something else... Best regards, Anthony Eyal on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer#post-5196 Sat, 08 Aug 2015 13:38:42 +0000 Eyal 5196@http://www.k-wave.org/forum/ Hi, I am very interested in this as well Please keep the updates coming :) Anthony on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer#post-5188 Fri, 31 Jul 2015 16:25:42 +0000 Anthony 5188@http://www.k-wave.org/forum/ It works in 3D (at least for a linear case, the nonlinear will come later...) :-) Have a good weekend, Anthony Anthony on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer#post-5187 Fri, 31 Jul 2015 14:43:05 +0000 Anthony 5187@http://www.k-wave.org/forum/ It appears to be fortuitous... When I use the velocity as source term in 2D, it does not match... I think the following observations are reliable (in 1D and 2D): - using only <code>p</code> as source term in the second simulation results in too low power. - using only <code>u</code> results in higher power (too much in 2D, the good one in 1D, probably by chance) - using both results in even higher power... Note : I always use Dirichlet source terms. However, I think I have a workaround: record the signal and set it as source term not only in one plane but in several consecutive planes (with pressure for example, I don't think it matters...). --> It works nicely for around 10 planes (source_width = 7 for example). I think it has a more "coercive effect" on the field in the second domain. I don't know if it makes sense considering the spatial colocation method, but it seems to work. I will try it in 3D and keep you informed. Best regards, Anthony <pre><code>clear all close all clc % Parameters CFL = 0.3; velocity_source = false; pressure_source = true; source_width = 7; % create the computational grid PMLsize = 20; Nx = 88; Ny = 88; dx = 0.1e-3; dy = 0.1e-3; kgrid = makeGrid(Nx, dx, Ny, dy); % define the properties of the propagation medium medium.sound_speed = 1500; medium.density = 1000; % create time array t_end = 1.1*sqrt(kgrid.x_size.^2 + kgrid.y_size.^2)/min(medium.sound_speed(:)); pixel_dim = min([kgrid.dx, kgrid.dy]); dt = CFL*pixel_dim/max(medium.sound_speed(:)); source_freq = 3e6; periode = 1/source_freq; NtimeStepsPerPeriod = ceil(periode/dt); dt = periode/NtimeStepsPerPeriod; kgrid.t_array = 0:dt:t_end; % create initial pressure distribution using makeDisc R = Ny*0.45; Rap = 28/38*R; source.p_mask = zeros(Nx, Ny); xc = round(Nx/2); yc = round(Ny/2); xgrid = (1:Nx)-xc; ygrid = ((1:Ny)-yc)'; rgrid = sqrt(bsxfun(@plus,xgrid.^2,ygrid.^2)); source.p_mask(logical(bsxfun(@times,(abs(rgrid-R)<1/2), bsxfun(@times,ygrid<0,abs(xgrid)<Rap)))) = 1; source.p_mode = 'dirichlet'; % define source source_mag = 0.1e6; source_cycles = 300; source_padding = 10; source.p = source_mag*sin(2*pi*kgrid.t_array*source_freq); % source.p = [source_mag * toneBurst(1/kgrid.dt, source_freq, source_cycles, 'Envelope', 'Rectangular'), zeros(1, source_padding)]; % define sensor sensor.mask = ones(Nx, Ny); sensor.record = {'p', 'u_non_staggered','u'}; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,... 'DisplayMask', 'off', 'PlotScale', [-0.5, 0.5]*source_mag,'PMLInside',false,'PlotPML',false,'PMLSize',[PMLsize PMLsize]); % Get the simulated pressure field over a plane planeIdx = 25; % arbitrary p_input = reshape(sensor_data.p,Nx,Ny,size(sensor_data.p,2)); p_input = p_input(planeIdx:planeIdx+source_width-1,:,:); p_input = reshape(p_input,[source_width*Ny,size(p_input,3)]); % Simulate the rest of the domain starting from the "intercepted plane" sourceShift = 10; Nx_new = Nx - planeIdx + 1 + sourceShift; kgrid2 = makeGrid(Nx_new, dx, Ny, dy); kgrid2.t_array = kgrid.t_array; % Set input if pressure_source source2.p_mask = zeros(Nx_new, Ny); source2.p_mask(1+sourceShift:1+sourceShift+source_width-1,:) = true; source2.p = p_input; source2.p_mode = 'dirichlet'; end if velocity_source ux_input = reshape(sensor_data.ux,Nx,Ny,size(sensor_data.ux,2)); ux_input = ux_input(planeIdx:planeIdx+source_width-1,:,:); ux_input = reshape(ux_input,[source_width*Ny,size(ux_input,3)]); uy_input = reshape(sensor_data.uy,Nx,Ny,size(sensor_data.uy,2)); uy_input = uy_input(planeIdx:planeIdx+source_width-1,:,:); uy_input = reshape(uy_input,[source_width*Ny,size(uy_input,3)]); if pressure_source source2.u_mask = source2.p_mask; else source2.u_mask = zeros(Nx_new, Ny); source2.u_mask(1+sourceShift:1+sourceShift+source_width-1,:) = true; end source2.ux = ux_input; source2.uy = uy_input; source2.u_mode = 'dirichlet'; end % set sensor2 sensor2.mask = ones(Nx_new, Ny); sensor2.record = {'p', 'u_non_staggered'}; % run the simulation sensor_data2 = kspaceFirstOrder2D(kgrid2, medium, source2, sensor2,... 'DisplayMask', 'off', 'PlotScale', [-0.5, 0.5]*source_mag,'PMLInside',false,'PlotPML',false,'PMLSize',[PMLsize PMLsize]); % time shift ux = timeShift(sensor_data.ux_non_staggered, kgrid.dt); uy = timeShift(sensor_data.uy_non_staggered, kgrid.dt); ux2 = timeShift(sensor_data2.ux_non_staggered, kgrid2.dt); uy2 = timeShift(sensor_data2.uy_non_staggered, kgrid2.dt); % Intensity NpointsToKeep = NtimeStepsPerPeriod; Ix = reshape(sensor_data.p.*ux,[Nx,Ny,size(ux,2)]); Iy = reshape(sensor_data.p.*uy,[Nx,Ny,size(uy,2)]); Ix = mean(Ix(:,:,end-NpointsToKeep+1:end),3); Iy = mean(Iy(:,:,end-NpointsToKeep+1:end),3); Ix2 = reshape(sensor_data2.p.*ux2,[Nx_new,Ny,size(ux2,2)]); Iy2 = reshape(sensor_data2.p.*uy2,[Nx_new,Ny,size(uy2,2)]); Ix2 = mean(Ix2(:,:,end-NpointsToKeep+1:end),3); Iy2 = mean(Iy2(:,:,end-NpointsToKeep+1:end),3); %% Display maskIdx = find(source.p_mask); [maskI,~] = ind2sub(size(source.p_mask),maskIdx); clear maskIdx endTrdI = max(maskI); figure(2) clf sidePowerLoss = ( cumsum(-Iy(:,1),1)+cumsum(Iy(:,end),1) ) * dx; compensatedPower = squeeze(sum(Ix,2)) * dy + sidePowerLoss; sidePowerLoss2 = ( cumsum(-Iy2(:,1),1) + cumsum(Iy2(:,end),1) ) * dx; compensatedPower2 = squeeze(sum(Ix2,2)) * dy + sidePowerLoss2 + sidePowerLoss(planeIdx); plot(compensatedPower, '-', 'LineWidth',2) hold on plot(planeIdx-sourceShift:planeIdx-sourceShift+numel(compensatedPower2)-1,compensatedPower2, '-', 'LineWidth',2) ylim([0 1.1*max([compensatedPower;compensatedPower2])]) ylabel('Integrated intensity [W]') xlabel('Main propagation axis [pixels]') set(gca,'FontSize',18) title(['Blue and red should be superimposed'])</code></pre> Anthony on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer#post-5186 Fri, 31 Jul 2015 10:49:20 +0000 Anthony 5186@http://www.k-wave.org/forum/ Hi, The daily additional comment has just arrived... and I think I got something :-) If you set the source term using (staggered) velocity instead of pressure, everything is perfect! I have no idea of the explanation but... at least, it seems to work! Anthony PS : by the way, it came that I looked at the code of the <code>timeShift()</code> function and I think you don't need <code>dt</code> as input parameter (it cancels in the computations). <pre><code>% Example to show recording and re-transmitting a signal using a Dirichlet % boundary condition % % author: Bradley Treeby % date: 19th July 2015 % last update: 28th July 2015 clear all close all clc source_type = 'u'; %'p' or 'u' or 'pu' % ========================================================================= % DEFINE LITERALS % ========================================================================= % perfectly matched layer pml_size = 20; pml_alpha = 2; % spatial grid Nx = 256 - 2*pml_size; dx = 0.1e-3; % temporal grid CFL = 0.3; t_end = 40e-6; % medium properties c0 = 1500; % source source_freq = 1.5e6; % [Hz] source_amp = 1e6; % [Pa] source_cycles = 50; source_padding = 100; % [time points] % positions source_pos = 1; % [grid points] split_pos = Nx/2; % [grid points] record_pos = Nx; % [grid points] % options plot_sim = false; % ========================================================================= % RUN FULL SIMULATION % ========================================================================= % create grid kgrid = makeGrid(Nx, dx); % define the medium properties medium.sound_speed = c0; medium.density = 1000; % define the time array t_array = makeTime(kgrid, c0, CFL, t_end); kgrid.t_array = t_array; % define the source mask source.p_mask = zeros(Nx, 1); source.p_mask(source_pos) = 1; % define the source source.p = [source_amp * toneBurst(1/kgrid.dt, source_freq, source_cycles, 'Envelope', 'Rectangular'), zeros(1, source_padding)]; source.p = filterTimeSeries(kgrid, medium, source.p); % define the sensor mask sensor.mask = ones(Nx, 1); % sensor.mask(record_pos, 1) = 1; % set the input arguments input_args = {'PMLInside', false, 'PMLSize', pml_size, 'PMLAlpha', pml_alpha,... 'PlotPML', false, 'PlotScale', [-1, 1]*source_amp, 'PlotSim', plot_sim}; % run the simulation sensor.record = {'p','u_non_staggered'}; full_sim_res = kspaceFirstOrder1D(kgrid, medium, source, sensor, input_args{:}); full_sim = full_sim_res.p(record_pos,:); % ========================================================================= % RUN SPLIT SIMULATION PART 1 % ========================================================================= % define sensor mask sensor.mask = zeros(Nx, 1); sensor.mask(split_pos, :) = 1; % run the simulation sensor.record = {'p','u'}; split_sim_part_1 = kspaceFirstOrder1D(kgrid, medium, source, sensor, input_args{:}); % ========================================================================= % RUN SPLIT SIMULATION PART 2 % ========================================================================= % take the previously recorded pressure and enforce as a Dirichlet boudnary % condition at the same position as it was recorded if strfind(source_type,'p') source.p_mask = zeros(Nx, 1); source.p_mask(split_pos, :) = 1; source.p = split_sim_part_1.p; source.p_mode = 'dirichlet'; end if strfind(source_type,'u') source.u_mask = zeros(Nx, 1); source.u_mask(split_pos, :) = 1; source.ux = split_sim_part_1.ux; source.u_mode = 'dirichlet'; end % define the sensor mask sensor.mask = ones(Nx, 1); % sensor.mask(record_pos, 1) = 1; sensor.record = {'p','u_non_staggered'}; % run the simulation split_sim_res = kspaceFirstOrder1D(kgrid, medium, source, sensor, input_args{:}); split_sim = split_sim_res.p(record_pos,:); % ========================================================================= %% DISPLAY % ========================================================================= % plot the signals figure(1); clf set(gcf,'Color',[1 1 1]) subplot(3, 1, 1); plot(full_sim * 1e-6, 'b-') hold on; plot(split_sim * 1e-6, 'r--') legend('Full Simulation', 'Split Simulation', 'Location', 'NorthWest'); title('Time Series'); xlabel('time index'); ylabel('pressure [MPa]'); % plot the difference subplot(3, 1, 2); plot(100 * abs(full_sim - split_sim) / max(full_sim(:))) title('Difference'); xlabel('time index'); ylabel('difference [% of max]'); % plot the power dt = kgrid.dt; full_power = full_sim_res.p.*timeShift(full_sim_res.ux_non_staggered,dt); split_power = split_sim_res.p.*timeShift(split_sim_res.ux_non_staggered,dt); full_power = mean(full_power(split_pos:end,end-99:end),2); split_power = mean(split_power(split_pos:end,end-99:end),2); subplot(3, 1, 3); set(gcf,'Color',[1 1 1]) plot(full_power,'LineWidth',2) hold on plot(split_power,'--','LineWidth',2) xlabel('spatial index (from split point)') ylabel('Intensity')</code></pre> Best regards, Anthony Anthony on "Simulation "layer by layer"" http://www.k-wave.org/forum/topic/simulation-layer-by-layer#post-5181 Thu, 30 Jul 2015 13:12:30 +0000 Anthony 5181@http://www.k-wave.org/forum/ Hi Brad, The amplitude and phase discrepancy seems to disappear for CFL = 1... I don't know what it means but maybe you will be able to interpret the clue :-) Anthony