function [hh,x,nn]=wellsim2d(q,L,m,K,Sy,dt,n,h)

// q is the discharge from well
// L is the length of lan under consideration
// for the well simulation, the area under consideration is 2D
//   and radial 'model' is used.
//   concentric circles/cells of annulus width 'dr' surround the well.
//   diameter of the largest radial is 'L'
// m is the number of cells
// K is conductance, Sy is the specific yield
// n is the number of time steps of simulation 
// dt is the time step for each iteration (round)
// h is the intial head of the cells measured at midpoints of cells
// hh: output, x:the X axis points

hh=h; 

// initialize number of iterations to 1, will go to nn (if required)
ii=1;     

pi=3.14;   // constant declaration

dr=L/m;    

//  cell(i) is the location (mid-point) of each cell, used for
// plotting purposes
for i=1:numCells
	// estimate the midpoint for each annulus
   cell(i)=dr/2+(i-1)*ringWidth;
end;


// start the iterations, till number of iterations done.
// each iteration, simulates head variation after time interval 'dt'

while ii<n

	// cell m is the well, cell 1 is the largest/farthest annulus
	// cell m and 1 are special cells, consider them later.
	for j=2:m-1

		// outflow from cell j to j+1
		// for the darcy's equation:
		// A = avg. head * circumference of annulus
	     
                head = (h(j+1)+h(j))/2; // head at inner surface of cell j
		innerR = ringWidth * (m-j);   // inner radius of cell j
		A = head*(2*pi*innerR); // A is the vertical inner surface area of the cylinder (cell j) 

		outflowRate = (h(j+1)-h(j))*K*A/ringWidth;

		// inflow from cell j-1 to cell j
                head = (h(j-1)+h(j))/2; // head at outer surface of cell j
		outerR = ringWidth * (m-j+1);  // outer radius of cell j 
		A = head*(2*pi*outerR); // A is the vertical outer surface area of the cylinder (cell j) 
		inflowRate = (h(j-1)-h(j))*K*A/ringWidth;

		netflowRate =  inflowRate + outflowRate;  
		netFlow =   netflowRate*dt;
		ringVolumeChange = netFlow/Sy;          

		r= (innerR+outerR)/2;		// average radius for annulus  		
		crossSectionArea = ringWidth*(2*pi*r); // Approximate horizontal cross section area for annulus
		dh = ringVolumeChange/crossSectionArea; 
		hn(j)=h(j)+dh;
	end;

	//cell m is the well, outflow from the well is 'q'
	head = (h(m)+h(m-1))/2; // head at outer surface of well
	r = ringWidth; // well radius
	A = head*pi*2*r;    // A is the vertical outer surface area of the well 

	inflowRate = (h(m-1)-h(m))*K*A/ringWidth;
	outflowRate = -q;
	netflowRate =  inflowRate + outflowRate;  
	netFlow =   netflowRate*dt;
	ringVolumeChange = netFlow/Sy;          

        crossSectionArea = (pi*r*r); // Well cross section is a complete circle, not an annulus, hence exact area = pi*r*r
	dh = ringVolumeChange/crossSectionArea;   

	hn(m) = h(m) + dh;

	// cell 1 is outermost annulus, no inflow only outflow
        head = (h(j+1)+h(j))/2; // head at inner surface of cell m
        innerR = ringWidth * (m-1);   // inner radius of cell 1
	A = head*(2*pi*innerR); // A is the vertical inner surface area of the cylinder (cell 1) 

	outflowRate = (h(2)-h(1))*K*A/ringWidth;

	inflowrate=0;

	netflowRate =  inflowRate + outflowRate;  
	netFlow =   netflowRate*dt;
	ringVolumeChange = netFlow/Sy;          

	r= (innerR+outerR)/2;		// average radius for annulus  		
	crossSectionArea = ringWidth*(2*pi*r); // Approximate horizontal cross section area for annulus
	dh = ringVolumeChange/crossSectionArea; 
	hn(m)=h(m)+dh;

	h=hn;
	ii=ii+1;

end;

hh=h;	// output heads of all annulus (note we are output a cross section)
nn=ii;  // output number of iterations

endfunction;

 

