Lap axisym

chunkie/+chnk/+axissymlap2d/g0funcall.m

@@ -0,0 +1,202 @@
+function [gvals, gdzs, gdrs, gdrps] = g0funcall(r, rp, dr, z, zp, dz, maxm)
+
+%
+% chnk.axissymlap2d.g0funcall evaluates a collection of axisymmetric Laplace
+% Green's functions, defined by the expression:
+%
+%     gfunc(n) = \frac{rp}{4\pi} * \int_0^{2\pi} 1/|x - x'| e^(-i n t) dt
+%
+% it is assumed that x = (x,0,z) otherwise the integral above should pick up a
+% phase factor out front of exp(i*n*phi), where phi is the azimuthal coordinate
+% of x in cylindrical coordinates.
+%
+% The extra factor of rp out front makes subsequent interfacing
+% with RCIP slightly easier. Modes 0 through maxm are returned, with gval(1) =
+% mode 0 and gval(maxm+1) = mode maxm. The function is even, so g_{-n} = g_n.
+%
+% The above scaling should be consistent with what is in
+% chnk.axissymlap2d.gfunc, which is for merely the zero-mode
+% 
+
+    twopi = 2*pi;
+    done = 1.0;
+
+    r0 = rp;
+    rzero = sqrt(r*r + r0*r0 + dz*dz);
+    alpha = 2*r*r0/rzero^2;
+    x = 1/alpha;
+    xminus = (dr*dr + dz*dz)/2/r/r0;
+
+    dxdr = (r^2 - r0^2 - (dz)^2)/2/r0/r^2;
+    dxdz = 2*(dz)/2/r/r0;
+    dxdr0 = (r0^2 - r^2 - (dz)^2)/2/r/r0^2;
+
+    %!
+    %! if xminus is very small, use the forward recurrence
+    %!
+
+    iffwd = 0;
+    if (x < 1.005)
+        iffwd = 1;
+        if ((x >= 1.0005d0) && (maxm > 163))
+            iffwd = 0;
+        end
+
+        if ((x >= 1.00005d0) && (maxm > 503))
+            iffwd = 0;
+        end
+
+        if ((x >= 1.000005d0) && (maxm > 1438))
+            iffwd = 0;
+        end
+
+        if ((x >= 1.0000005d0) && (maxm > 4380))
+            iffwd = 0;
+        end
+
+        if ((x >= 1.00000005d0) && (maxm > 12307))
+            iffwd = 0;
+        end
+    end    
+
+    gvals = zeros(maxm+1,1);
+    gdzs = zeros(maxm+1,1);
+    gdrs = zeros(maxm+1,1);
+    gdrps = zeros(maxm+1,1);
+
+    if (iffwd == 1)
+
+        [q0, q1, dq0] = chnk.axissymlap2d.qleg_half(xminus);
+        dq1 = (-q0 + x*q1)/2/(x+1)/xminus;
+
+        half = done/2;
+
+        fac = sqrt(rp/r)/twopi;
+        gvals(1) = fac*q0;
+        gvals(2) = fac*q1;
+
+        derprev = fac*dq0;
+        der = fac*dq1;
+
+        % the z derivatives
+        gdzs(1) = derprev*dxdz;
+        gdzs(2) = der*dxdz;
+
+        % the r derivatives
+        gdrs(1) = fac*(dq0*dxdr - q0/2/r);
+        gdrs(2) = fac*(dq1*dxdr - q1/2/r);
+
+        % the rp derivatives
+        gdrps(1) = fac*(dq0*dxdr0 - q0/2/r0);
+        gdrps(2) = fac*(dq1*dxdr0 - q1/2/r0);
+
+        % run upward recursion for the Q's to calculate them things
+        for i = 1:(maxm-1)
+            gvals(i+2) = (2*i*x*gvals(i+1) - (i-half)*gvals(i))/(i+half);
+            dernext = (2*i*(gvals(i+1)+x*der) - (i-half)*derprev)/(i+half);
+            gdrs(i+2) = (dernext*dxdr - gvals(i+2)/2/r);
+            gdzs(i+2) = dernext*dxdz;
+            gdrps(i+2) = (dernext*dxdr0 - gvals(i+2)/2/r0);
+            derprev = der;
+            der = dernext;
+        end
+
+        return
+    end
+
+
+    %!
+    %! if here, xminus > .005, so run forward and backward recurrence
+    %!
+
+    %!
+    %! run the recurrence, starting from maxm, until it has exploded
+    %! for BOTH the values and derivatives
+    %!
+    done = 1;
+    half = done/2;
+    f = 1;
+    fprev = 0;
+    der = 1;
+    derprev = 0;
+    maxiter = 100000;
+    upbound = 1.0e19;
+
+    for i = maxm:maxiter
+        fnext = (2*i*x*f - (i-half)*fprev)/(i+half);
+        dernext = (2*i*(x*der+f) - (i-half)*derprev)/(i+half);
+        if (abs(fnext) >= upbound) 
+            if (abs(dernext) >= upbound)
+                nterms = i+1;
+                break
+            end
+        end
+        fprev = f;
+        f = fnext;
+        derprev = der;
+        der = dernext;
+    end
+
+    %!
+    %! now start at nterms and recurse down to maxm
+    %!
+    if (nterms < 10)
+        nterms = 10;
+    end
+
+    fnext = 0;
+    f = 1;
+    dernext = 0;
+    der = 1;
+
+    % run the downward recurrence
+    for j = 1:(nterms-maxm+1)
+        i = nterms-j+1;
+        fprev = (2*i*x*f - (i+half)*fnext)/(i-half);
+        fnext = f;
+        f = fprev;
+        derprev = (2*i*(x*der+f) - (i+half)*dernext)/(i-half);
+        dernext = der;
+        der = derprev;
+    end
+
+    gvals(maxm) = f;
+    gvals(maxm+1) = fnext;
+
+    ders = zeros(maxm+1,1);
+    ders(maxm) = der;
+    ders(maxm+1) = dernext;
+
+    for j = 1:(maxm-1)
+        i = maxm-1-j+1;
+        gvals(i) = (2*i*x*gvals(i+1) - (i+half)*gvals(i+2))/(i-half);
+        ders(i) = (2*i*(x*ders(i+1)+gvals(i+1)) - (i+half)*ders(i+2))/(i-half);
+   end
+
+
+    % !
+    % ! normalize the values, and use a formula for the derivatives
+    % !
+    [q0, q1, dq0] = chnk.axissymlap2d.qleg_half(xminus);
+    dq1 = (-q0 + x*q1)/2/(x+1)/xminus;
+
+    ratio = q0/gvals(1)*sqrt(rp/r)/twopi;
+
+    for i = 1:(maxm+1)
+        gvals(i) = gvals(i)*ratio;
+    end
+
+    ders(1) = dq0*sqrt(rp/r)/twopi;
+    ders(2) = dq1*sqrt(rp/r)/twopi;
+    for i = 2:maxm
+       ders(i+1) = -(i-.5d0)*(gvals(i) - x*gvals(i+1))/(1+x)/xminus;
+    end
+
+    %
+    % and scale the gradients properly everyone...
+    %
+    gdzs = ders*dxdz;
+    gdrs = ders*dxdr - gvals/2/r;
+    gdrps = ders*dxdr0 - gvals/2/r0;
+
+end

chunkie/+chnk/+axissymlap2d/g0funcall_vec.m

@@ -0,0 +1,141 @@
+function [gval, gdz, gdr, gdrp, gdrpr, gdzz, gdrz, gdrpz] = g0funcall_vec(r, rp, dr, z, zp, dz, m)
+%
+% chnk.axissymlap2d.g0funcall evaluates a collection of axisymmetric Laplace
+% Green's functions, defined by the expression:
+%
+%     gfunc(n) = pi*rp * \int_0^{2\pi} 1/|x - x'| e^(-i n t) dt
+%
+% Modes 0 through maxm are returned, with gval(1) = mode 0 and
+% gval(maxm+1) = mode maxm.
+%
+
+    r  = reshape(r, [1,size(r,1),size(r,2)]);
+    rp = reshape(rp,[1,size(rp,1),size(rp,2)]);
+    dr = reshape(dr,[1,size(dr,1),size(dr,2)]);
+    z  = reshape(z, [1,size(z,1),size(z,2)]);
+    zp = reshape(zp,[1,size(zp,1),size(zp,2)]);
+    dz = reshape(dz,[1,size(dz,1),size(dz,2)]);
+
+    % we handle r=0 and rp=0 cases separetly 
+    src_on_axis = (rp == 0);
+    targ_on_axis = (r == 0) & (rp ~= 0);
+
+    r_safe = r;
+    rp_safe = rp;
+
+    % dummy values whos results will be overwritten
+    r_safe(r_safe == 0) = 1;
+    rp_safe(rp_safe == 0) = 1;
+
+    % case: rp != 0 and r != 0
+    t = (dz.^2 + dr.^2)./(2.*r_safe.*rp_safe);
+    chi = t + 1;
+
+    [qm, qmd, qmdd] = chnk.axissymlap2d.qleg_half_miller_vec(t,m);
+
+    prefac = sqrt(rp_safe./r_safe)/(2*pi);
+
+    gval = prefac.*qm;
+
+    gdz  = prefac.*qmd ...
+          ./(rp_safe.*r_safe).*dz;
+
+    rfac = -r_safe/2.*qm + (-(1+t).*r_safe + rp_safe).*qmd;
+    gdrp = prefac./(rp_safe.*r_safe).*rfac;
+
+    rfac = -rp_safe/2.*qm + (-(1+t).*rp_safe + r_safe).*qmd;
+    gdr  = prefac./(rp_safe.*r_safe).*rfac;
+
+    rfac = 1./(rp_safe.*r_safe).*qmd ...
+        + (dz./(rp_safe.*r_safe)).^2.*qmdd;
+    gdzz = prefac.*rfac;
+
+    rfac = -3./(2*r_safe.^2.*rp_safe).*qmd ...
+        + (-chi./(r_safe.^2.*rp_safe) + 1./(r_safe.*rp_safe.^2)).*qmdd;
+    gdrz = prefac.*dz.*rfac;
+
+    rfac = -3./(2*rp_safe.^2.*r_safe).*qmd ...
+        + (-chi./(rp_safe.^2.*r_safe) + 1./(rp_safe.*r_safe.^2)).*qmdd;
+    gdrpz = prefac.*dz.*rfac;
+
+    rfac = 1./(4*r_safe.*rp_safe).*qm ...
+        + (2*chi./(rp_safe.*r_safe) ...
+        - 3./(2*r_safe.^2) ...
+        - 3./(2*rp_safe.^2)).*qmd ...
+        + (-chi./r_safe + 1./rp_safe).*(-chi./rp_safe + 1./r_safe).*qmdd;
+    gdrpr = prefac.*rfac;
+
+    % case: rp = 0
+    if any(src_on_axis(:))
+        gval(:,:,src_on_axis) = 0;
+        gdz(:,:,src_on_axis) = 0;
+        gdr(:,:,src_on_axis) = 0;
+        gdrp(:,:,src_on_axis) = 0;
+        gdrpr(:,:,src_on_axis) = 0;
+        gdzz(:,:,src_on_axis) = 0;
+        gdrz(:,:,src_on_axis) = 0;
+        gdrpz(:,:,src_on_axis) = 0;
+    end
+
+    % case: r = 0
+    if any(targ_on_axis(:))
+        gval(:,targ_on_axis) = 0;
+        gdz(:,targ_on_axis) = 0;
+        gdr(:,targ_on_axis) = 0;
+        gdrp(:,targ_on_axis) = 0;
+        gdrpr(:,targ_on_axis) = 0;
+        gdzz(:,targ_on_axis) = 0;
+        gdrz(:,targ_on_axis) = 0;
+        gdrpz(:,targ_on_axis) = 0;
+
+        rp0 = rp(targ_on_axis);
+        dz0 = dz(targ_on_axis);
+
+        tmp = gval(1,:,:);
+        tmp(targ_on_axis) = rp0 ...
+            ./ sqrt(rp0.^2 + dz0.^2)/2;
+        gval(1,:,:) = tmp;
+
+        tmp = gdz(1,:,:);
+        tmp(targ_on_axis) = -rp0.*dz0 ...
+            ./ sqrt(rp0.^2 + dz0.^2).^3/2;
+        gdz(1,:,:) = tmp;
+
+        tmp = gdrp(1,:,:);
+        tmp(targ_on_axis) = -rp0.^2 ...
+            ./ sqrt(rp0.^2 + dz0.^2).^3/2;
+        gdrp(1,:,:) = tmp;
+
+        tmp = gdzz(1,:,:);
+        tmp(targ_on_axis) = rp0.* ( ...
+                -1 ./ sqrt(rp0.^2 + dz0.^2).^3 ...
+                + 3*dz0.^2 ...
+                ./ sqrt(rp0.^2 + dz0.^2).^5 ...
+            )/2;
+        gdzz(1,:,:) = tmp;
+
+        tmp = gdrpz(1,:,:);
+        tmp(targ_on_axis) = 3*rp0.^2.*dz0 ...
+            ./ sqrt(rp0.^2 + dz0.^2).^5/2;
+        gdrpz(1,:,:) = tmp;
+
+        if m >= 1
+            tmp = gdr(2,:,:);
+            tmp(targ_on_axis) = rp0.^2 ...
+                ./ sqrt(rp0.^2 + dz0.^2).^3/4;
+            gdr(2,:,:) = tmp;
+
+            tmp = gdrz(2,:,:);
+            tmp(targ_on_axis) = -3*rp0.^2.*dz0 ...
+                ./ sqrt(rp0.^2 + dz0.^2).^5/4;
+            gdrz(2,:,:) = tmp;
+
+            tmp = gdrpr(2,:,:);
+            tmp(targ_on_axis) = ( ...
+                    rp0 ./ sqrt(rp0.^2 + dz0.^2).^3 ...
+                    - 3*rp0.^3 ./ sqrt(rp0.^2 + dz0.^2).^5 ...
+                )/4;
+            gdrpr(2,:,:) = tmp;
+        end
+    end
+end

chunkie/+chnk/+axissymlap2d/gaus_agm.m

@@ -0,0 +1,38 @@
+function [rk0,re0] = gaus_agm(x)
+
+eps  = 1E-15;
+a    = sqrt(2./(x+1));
+delt = 1./sqrt(1-a.*a);
+aa0  = delt + sqrt(delt.*delt-1);
+bb0   = 1./(delt+sqrt(delt.*delt-1));
+a0  = ones(size(delt));
+b0  = 1./delt;
+
+
+fact = ((a0+b0)/2).^2;
+
+for i=1:1000
+    a1 = (a0+b0)/2;
+    b1 = sqrt(a0.*b0);
+
+    aa1 = (aa0+bb0)/2;
+    bb1 = sqrt(aa0.*bb0);
+    a0 = a1;
+    b0 = b1;
+    aa0 = aa1;
+    bb0 = bb1;
+
+    c0 = (a1-b1)/2;
+    fact = fact-(c0.*c0)*2^(i);
+    drel = abs(a0-b0)./abs(a0);
+    drel2= abs(aa0-bb0)./abs(aa0);
+    if (max(drel+drel2) <2*eps)
+        break
+    end
+
+end
+
+rk0 = pi./(2*aa0.*sqrt(1-a.*a));
+re0 = pi*fact./(2*a0);
+
+end