function [phi2,DD] = sweep3(phi2,n,h,N,phi,iix,iiy) %% second order correction %% stencil is chosen using the determinant %% the bigger the det. the more upwind ??? %% loop on grid point according to sweeping DD = zeros(size(phi)) ; for j = iiy for i = iix %%%% shutof source data (exact) if (phi(i,j) > 0.1) %% find the second order stencil using first order solution %% upwind quadrant does not work for grid lines %% if ( phi(i-1, j) < phi(i,j) ) ix = -1 ; a1 = phi(i,j)-phi(i-1,j) ; elseif ( phi(i+1, j) < phi(i,j) ) ix = 1 ; a1 = phi(i,j)-phi(i+1,j) ; else ix = 0 ; a1 = 0 ; end if ( phi(i, j-1) < phi(i,j) ) b1 = phi(i,j) - phi(i,j-1) ; jy = -1 ; elseif ( phi(i,j) > phi(i,j+1) ) b1 = phi(i,j) - phi(i,j+1) ; jy = 1 ; else b1 = 0 ; jy = 0 ; end alpha = max(eps,b1/max(eps,a1)) ; %% this is the tangent parameter (abs().. always positive %%% central stencil fx = phi2(i+ix,j) ; fy = phi2(i,j+jy) ; fxy = phi2(i+ix,j+jy) ; %% right X upwind d2 = phi2(i+ix,j+1)+phi2(i+ix,j-1)-2*fx ; B = 1/2*(phi2(i+ix,j+1)-phi2(i+ix,j-1)) + jy*alpha*d2 ; A = -fx + 1/2*(alpha)^2*d2 ; c1r = 1 ; c2r = 2*A ; c3r = B^2 + A^2 - h^2*n(i,j)^2 ; Dr = c2r^2-4*c1r*c3r ; %% left Y upwind d2 = phi2(i+1,j+jy)+phi2(i-1,j+jy)-2*fy ; A = 1/2*(phi2(i+1,j+jy)-phi2(i-1,j+jy)) + ix/alpha*d2 ; B = -fy + 1/2*(1/alpha)^2*d2 ; c1l = 1 ; c2l = 2*B ; c3l = B^2 + A^2 - h^2*n(i,j)^2 ; Dl = c2l^2-4*c1l*c3l ; %% center A = (alpha/2*(fy-fxy) - (1-alpha/2)*fx) ; B = (1/(2*alpha)*(fx-fxy) - (1-1/(2*alpha))*fy) ; c1c = (1-alpha/2)^2 + (1-1/(2*alpha))^2 ; c2c = 2*A*(1-alpha/2) + 2*B*(1-1/(2*alpha)) ; c3c = B^2 + A^2 - h^2*n(i,j)^2 ; Dc = c2c^2-4*c1c*c3c ; %% stencil selected on the max Optimization needed ! D = max([Dc Dr Dl]) ; DD(i,j) = D ; if ( D >=0) if ( Dc == D) phi2(i,j) =(-c2c+sqrt(Dc))/(2*c1c) ; elseif (Dl == D) phi2(i,j) =(-c2l+sqrt(Dl))/(2*c1l) ; else %% c'est Dr phi2(i,j) =(-c2r+sqrt(Dr))/(2*c1r) ; end end end %% for end %% for end %% if