spherical_lap

 calculates laplacian on spherical surface

0001 function v_lap = spherical_lap(welec,x,y,z,xs,ys,zs,p,q)
0002 %
0003 % calculates laplacian on spherical surface
0004 %
0005 welec = welec.^2;
0006 n = size(x,2);
0007 ns = size(xs,2);
0008 [azs,els,rs] = cart2sph(xs,ys,zs);
0009 els = (pi/2)*ones(size(els,1),size(els,2)) - els;
0010 [az,el,r] = cart2sph(x,y,z);
0011 el = (pi/2)*ones(size(el,1),size(el,2)) - el;
0012 
0013 % calculate trig functions to keep things easier
0014 
0015 st = sin(els);
0016 ct = cos(els);
0017 sp = sin(azs);
0018 cp = cos(azs);
0019 s2t = sin(2*els);
0020 c2t = cos(2*els);
0021 s2p = sin(2*azs);
0022 c2p = cos(2*azs);
0023 
0024 % calculate utta uppa utt upp ucc
0025 
0026 % the real enchilada is uuxyz
0027 %
0028 % note for future:  this is where i have to rederive for arbitrary surface normals
0029 %
0030 
0031 uuxyz = 2*q(4)+2*q(6)+2*q(10) - (2*st.*(q(2)*cp+q(3)*sp)./rs +2*q(7)*ct./rs+6*(st.^2).*(q(4)*(cp.^2)+q(6)*(sp.^2)+q(5)*sp.*cp)+6*st.*ct.*(q(8)*cp+q(9)*sp)+6*q(10)*ct.^2);
0032 
0033 
0034 %
0035 % its the unfortunate time where we deal witht logarithmic mess
0036 %
0037 
0038 %
0039 % zero registers
0040 %
0041 
0042 smpp = zeros(size(st,1),size(st,2));
0043 smtt = zeros(size(st,1),size(st,2));
0044 smp = zeros(size(st,1),size(st,2));
0045 smt = zeros(size(st,1),size(st,2));
0046 ttcomp = zeros(size(st,1),size(st,2));
0047 rrcomp = zeros(size(st,1),size(st,2));
0048 
0049 for j = 1:n
0050      a = r(j)*(st.*cp-sin(el(j))*cos(az(j))*ones(size(st,1),size(st,2)));
0051     b = r(j)*(st.*sp-sin(el(j))*sin(az(j))*ones(size(st,1),size(st,2)));
0052     c = r(j)*(ct-cos(el(j))*ones(size(st,1),size(st,2)));
0053     str = a.^2+b.^2+c.^2;
0054     strw = str+welec*ones(size(st,1),size(st,2));
0055     comterm = 4*str./strw-((str./strw).^2)+2*log(strw);
0056     comterm2 = 2*(2*str.*log(strw)+(str.^2)./strw);
0057     tcomp = 3*comterm2+4*str.*comterm;
0058     dr = 2*(a.*st.*cp+b.*st.*sp+c.*ct);
0059     d2r2 = 2;
0060     rcomp = dr.*comterm2+d2r2*r(j)*comterm2/2+r(j)*(dr.^2).*comterm;
0061     ttcomp = ttcomp + p(j)*tcomp;
0062     rrcomp = rrcomp + p(j)*rcomp/r(j);
0063 end;
0064 v_lap = - (ttcomp+uuxyz-rrcomp);
0065 
0066 
0067 
0068 
0069 
0070 
0071 
0072 
0073 
0074