% Computes the potential around a three-phase high tension line, 20 m
% above ground, with 5 m horizontal separation between the three lines. 
% Artificial boundary at x = +-30 m and 50 m above ground.  Uses the
% method described in Sect. 12.1 in the textbook by N.N Rao, and 
% computes only one iteration of the scheme. 
% 
% Input to this script:
% Dx:   no. sampling points between conductors
% iter: iteration counter, has to be set to zero initially
% istp: number of iterations between each plot
%
if ~exist('Dx');      Dx = 2;  end
if ~exist('iter');  iter = 0;  end
if ~exist('istp');  istp = 1;  end
if iter == 0     % Initialization
   dh =  5./Dx;  % distance between sampling points in m
   Kx =  6.*Dx;
   x  =  dh.*(-Kx:Kx);                 % Set up x-axis
   Lx =  2*Kx;
   Ly =  10.*Dx;
   y  =  dh.*(0:Ly);                   % Set up y-axis
   Dy =  4.*Dx;                        % distance above ground
   V0 =  [-100000 0 100000].*sqrt(3);  % Tension on the three phases
%  V0 =  [-100000 200000 -100000];     % Quarter-period delayed tension 
   fi =  zeros(Ly+1,Lx+1);             % Initial potential distribution
   fi(Dy+1,Kx+1-Dx) =  V0(1); 
   fi(Dy+1,Kx+1)    =  V0(2); 
   fi(Dy+1,Kx+1+Dx) =  V0(3); 
end
iter  =  iter + istp;
for iters = 1:istp
    fi(2:Ly,2:Lx) = 0.25 .* ...
    (fi(1:(Ly-1),2:Lx) + fi(3:(Ly+1),2:Lx) + ...
    fi(2:Ly,1:(Lx-1)) + fi(2:Ly,3:(Lx+1)));
    fi(:,1)     =  fi(:,2);         % put normal derivatives to zero
    fi(:,Lx+1)  =  fi(:,Lx);
    fi(Ly+1,:)  =  fi(Ly,:);
    fi(Dy+1,Kx-Dx+1) =  V0(1);      % reset line voltages
    fi(Dy+1,Kx+1)    =  V0(2); 
    fi(Dy+1,Kx+Dx+1) =  V0(3); 
end
fiplot   =  max(fi,-20000);    
fiplot   =  min(fiplot,20000);  % avoid plotting very high voltages
contour (x,y,fiplot,20);


