#include <math.h>
#include <stdio.h>
#include <Riostream.h>

#include "TMath.h"
#include "TCanvas.h"
#include "TMatrixD.h" 

#define  COLUMNS       201     // Must be ODD, Must be Multiple of 200
#define  ROWS         2001     // Must be ODD
#define  WIREX1          0
#define  WIREX2          0
#define  WIREX3          0
#define  WIREY1        200
#define  WIREY2        400
#define  WIREY3       1000
#define  V0              0     // Pad Readout Plane (grounded)
#define  V1          -1170     // Anode potential (should be positive)
#define  V2              0     // Cathode Grid (grounded)
#define  V3             40     // Gated Grid I  (-40 or -115 or -190)
#define  V4            190     // Gated Grid II (-190 or -115 or -40)
#define  V8            246     // -380Volts on pseudo CM at 3 cm
#define  RADII0          8     // Square wires (for now)
#define  RADII1          4     // Square wires (for now)
#define  ITERATIONS   1000

void Many( ) 
{

  TCanvas *c1, *c2, *c3;        // Canvases for drawing graphs
  c1 = new TCanvas("c1", "Voltage Plot", 100, 100, 800, 500) ;
  c2 = new TCanvas("c2","Ex Field Plot", 150, 150, 800, 500) ;
  c3 = new TCanvas("c3","Ey Field Plot", 200, 200, 800, 500) ;

  TMatrixD BigV(ROWS,COLUMNS), BigEx(ROWS-1,COLUMNS-1), BigEy(ROWS-1,COLUMNS-1) ;

  Float_t t_nagle   = TMath::Cos(TMath::Pi()/ROWS) + TMath::Cos(TMath::Pi()/COLUMNS) ;  // Nagles's overrelaxation parameter
  Float_t OVERRELAX = ( 8 - TMath::Sqrt(64-16*t_nagle*t_nagle)) / (t_nagle*t_nagle)  ;  // Nagles's overrelaxation parameter

  for ( Int_t i = 0 ; i < ROWS ; i++ )
    {
      for ( Int_t j = 0 ; j < COLUMNS ; j++ )
	{

	  BigV(i,j) =  V8 * (float)i/(ROWS-1)                                        ;  // Fill array with approximate solution
	  if ( i==0 )  BigV(i,j) = V0                                                ;  // Boundary Conditions
	  if ( i==(ROWS-1) )  BigV(i,j) = V8                                         ;  // Boundary Conditions
	  Int_t m = j%400 ; if ( m > 200 ) m = m - 400 ;                             ;  // help select a wire every 4 millimeters (grid specific)
	  Int_t n = j%100 ; if ( n > 50  ) n = n - 100 ;                             ;  // help select a wire every 1 millimeter  (grid specific)
  	  if ( (i>=(WIREY1-RADII0)) && (i<=(WIREY1+RADII0)) && (m>=(WIREX1-RADII0)) && (m<=(WIREX1+RADII0)) ) BigV(i,j)=V1 ; // Special condtions for wires
	  if ( (i>=(WIREY2-RADII1)) && (i<=(WIREY2+RADII1)) && (n>=(WIREX2-RADII1)) && (n<=(WIREX2+RADII1)) ) BigV(i,j)=V2 ; // Special condtions for wires
	  if ( (i>=(WIREY3-RADII1)) && (i<=(WIREY3+RADII1)) && (n>=(WIREX3-RADII1)) && (n<=(WIREX3+RADII1)) )
	    {
	      if ( ((j+50)/100)%2 == 0 ) BigV(i,j) = V3   ;
	      else BigV(i,j) = V4                    ;
	    }
	}
    }

  for ( Int_t k = 0 ; k < ITERATIONS ; k++ )
    {
      for ( Int_t i = 1 ; i < (ROWS-1) ; i++ )
	{
	  for ( Int_t j = 0 ; j < COLUMNS ; j++ )
	    {
	      Int_t m = j%400 ; if ( m > 200 ) m = m - 400 ;                            ;  // help select a wire every 4 millimeters (grid specific)
	      Int_t n = j%100 ; if ( n > 50  ) n = n - 100 ;                            ;  // help select a wire every 1 millimeter  (grid specific)
	      /*	      if ( (i>=(WIREY1-RADII0)) && (i<=(WIREY1+RADII0)) && (m>=(WIREX1-RADII0)) && (m<=(WIREX1+RADII0)) )
		{
		  if      ( i==(WIREY1+RADII0) && m==(WIREX1+RADII0) ) ;
		  else if ( i==(WIREY1-RADII0) && m==(WIREX1+RADII0) ) ;
		  else if ( i==(WIREY1+RADII0) && m==(WIREX1-RADII0) ) ;
		  else if ( i==(WIREY1-RADII0) && m==(WIREX1-RADII0) ) ;
		  else    continue ; // Special condtions for wires
		  }*/
	      
	      if ( (i>=(WIREY1-RADII0)) && (i<=(WIREY1+RADII0)) && (m>=(WIREX1-RADII0)) && (m<=(WIREX1+RADII0)) ) continue ; // Special conditions for wires
	      if ( (i>=(WIREY2-RADII1)) && (i<=(WIREY2+RADII1)) && (n>=(WIREX2-RADII1)) && (n<=(WIREX2+RADII1)) ) continue ; // Special conditions for wires
	      if ( (i>=(WIREY3-RADII1)) && (i<=(WIREY3+RADII1)) && (n>=(WIREX3-RADII1)) && (n<=(WIREX3+RADII1)) ) continue ; // Special conditions for wires
	      if ( j == 0 )                                                             // Solve Laplace's equation with a Reflection Boundary Condition
		BigV(i,0) = OVERRELAX * 0.25 * ( BigV(i-1,0) + BigV(i,1) + BigV(i+1,0) + BigV(i,1) ) - (OVERRELAX-1)*BigV(i,0) ; 
	      else if ( j == (COLUMNS-1) )                                              // Solve Laplace's equation with a Reflection Boundary Condition
		BigV(i,COLUMNS-1) = OVERRELAX * 0.25 * ( BigV(i-1,COLUMNS-1) + BigV(i,COLUMNS-2) + BigV(i+1,COLUMNS-1) + BigV(i,COLUMNS-2) ) 
		  - (OVERRELAX-1)*BigV(i,COLUMNS-1) ;          
	      else                                                                      // Solve LaPlace's equation at an interior point
		BigV(i,j) = OVERRELAX * 0.25 * ( BigV(i-1,j) + BigV(i,j+1) + BigV(i+1,j) + BigV(i,j-1) ) - (OVERRELAX-1)*BigV(i,j)   ;
	    }
	}
    }

  for ( Int_t i = 0 ; i < (ROWS-1) ; i++ )
    {
      for ( Int_t j = 0 ; j < (COLUMNS-1) ; j++ )
	{
	  BigEx(i,j) = 1000 * 0.5 * ( BigV(i,j+1) - BigV(i,j) + BigV(i+1,j+1) - BigV(i+1,j) ) ;
	  BigEy(i,j) = 1000 * 0.5 * ( BigV(i+1,j) - BigV(i,j) + BigV(i+1,j+1) - BigV(i,j+1) ) ;
	}
    }

  c1 -> cd() ;
  BigV.Draw("surf2") ;
  c1 -> Update() ;
  c2 -> cd() ;
  BigEx.Draw("surf2") ;
  c2 -> Update() ;
  c3 -> cd() ;
  BigEy.Draw("surf2") ;
  c3 -> Update() ;

}