clear; lw=1;
a=2; b=1.4; c=sqrt(a^2-b^2); n=200;
theta=linspace(0,2*pi,n+1);
theta=[theta(end-1),theta];
x=a*cos(theta);
y=b*sin(theta);
ep=0.1;
A1=[-1.2;0]; A2=[1.2;0];
p=[x;y];
cal_Nvec = @(pp,the) [cos(the),-sin(the);sin(the),cos(the)]*pp;
temp_p1=p(:,3:end)-p(:,2:end-1);
temp_p1=temp_p1./vecnorm(temp_p1);
Np1=cal_Nvec(temp_p1,pi/2);
temp_p2=p(:,1:end-2)-p(:,2:end-1);
temp_p2=temp_p2./vecnorm(temp_p2);
Np2=cal_Nvec(temp_p2,-pi/2);
dir_p=(Np1+Np2)/2;
dir_p=dir_p./vecnorm(dir_p);
f=1./(abs(vecnorm(A1-p(:,2:end-1))-vecnorm(A2-p(:,2:end-1)))+ep);
fdir=f.*dir_p;
figure(1); clf; hold on;
plot(p(1,:),p(2,:),'ko-')
plot(A1(1),A1(2),'ro',A2(1),A2(2),'ro')
quiver(p(1,2:end-1),p(2,2:end-1),fdir(1,:),fdir(2,:));
axis image
axis off

 

 

 

 

 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 

clear; lw=1;
a=10; b=7; c=7; n=20;
theta=linspace(0,2*pi,n+1);
x=zeros(n,n+1);
y=x; z=x;
z(:,1)=linspace(-c,c,n)';
z(1,:)=z(1,1); z(n,:)=z(n,1);
for k=2:n-1
x(k,:)=a*sqrt(1-z(k)^2/c^2)*cos(theta);
y(k,:)=b*sqrt(1-z(k)^2/c^2)*sin(theta);
z(k,:)=z(k,1);
end
x=[x(:,end-1),x];
x=[x(2,:);x;x(end-1,:)];
y=[y(:,end-1),y];
y=[y(2,:);y;y(end-1,:)];
z=[z(:,end-1),z];
z=[z(2,:);z;z(end-1,:)];
p=zeros(3,n+2,n+2);
p(1,:,:)=x;
p(2,:,:)=y;
p(3,:,:)=z;
Nt=2;
ep=0.1;
A1=[-6;0;0]; A2=[6;0;0];
%%% 각 p의 방향벡터와 크기 결정 %%%
[Nx,Ny,Nz] = surfnorm(x(2:end-1,2:end-1),y(2:end-1,2:end-1),z(2:end-1,2:end-1));
Nx=-Nx;
Ny=-Ny;
Nz=-Nz;
Nxx=Nx./sqrt(Nx.^2+Ny.^2+Nz.^2+eps);
Nyy=Ny./sqrt(Nx.^2+Ny.^2+Nz.^2+eps);
Nzz=Nz./sqrt(Nx.^2+Ny.^2+Nz.^2+eps);
f=zeros(n);
for i=1:n
for j=1:n
f(i,j)=1./(abs(vecnorm(A1-p(:,i+1,j+1))-vecnorm(A2-p(:,i+1,j+1)))+ep);
fdir(:,i,j)=[Nxx(i,j);Nyy(i,j);Nzz(i,j)]*f(i,j);
end
end
%%% 결과 그리기 %%%
figure(1); clf; hold on; view(10,10);
%surf(x(2:end,2:end),y(2:end,2:end),z(2:end,2:end),'EdgeColor','none','FaceColor','k','FaceAlpha',0.01)
plot3(A1(1),A1(2),A1(3),'ro',A2(1),A2(2),A2(3),'ro')
camlight('left');
quiver3(p(1,2:end-1,2:end-1),p(2,2:end-1,2:end-1),p(3,2:end-1,2:end-1),fdir(1,:,:),fdir(2,:,:),fdir(3,:,:));
axis image; axis off;