// TEST SCRIPT FOR GetContourList Function
// Author: Josh de Bever
//         CSI Medical Physics Group
//         The University of Western Ontario
//         London, Ontario, Canada
//   Date: Aug. 26, 2004


CERNtest(){
    gROOT->Reset();
    gSystem->Load("libphysics");
 
    const Double_t PI = TMath::Pi(); 
    
    TCanvas* c = new TCanvas("c","Sample",0,0,700,600);
 
    Int_t i, j, pti, TotalConts;
    
    Int_t nZsamples = 80;
    Int_t nPhiSamples = 80;

    Double_t HofZwavelength = 4.0; // 4meters 
    Double_t dZ =  HofZwavelength/(Double_t)(nZsamples - 1);
    Double_t dPhi = 2*PI/(Double_t)(nPhiSamples - 1);
 
    TArrayD z(nZsamples);
    TArrayD HofZ(nZsamples);
    TArrayD phi(nPhiSamples);
    TArrayD FofPhi(nPhiSamples);

    TGraph2D *streamFn = new TGraph2D();
  
    // Discritize Z and Phi Values
    for ( i = 0; i < nZsamples; i++) {
	z[i] = (i)*dZ - HofZwavelength/2.0;
	
	HofZ[i] = SawTooth(z[i], HofZwavelength);
	printf("%d: Z = % f  H(Z) = % f\n", i, z[i], HofZ[i]);
    }

    for(Int_t i=0; i < nPhiSamples; i++){
	phi[i] = (i)*dPhi;      
	FofPhi[i] = sin(phi[i]); 
    }

    printf("Input Functions Calculated\n");
     
    // Load Graph2D
    for (Int_t i = 0; i < nZsamples; i++) {
	for(Int_t j = 0; j < nPhiSamples; j++){

	    streamFn->SetPoint(pti , z[i], phi[j], HofZ[i]*FofPhi[j]);
	    pti++;

	} 
    }

    gStyle->SetPalette(1); 
    streamFn->Draw("surf1");

    // *************************************************************************
    // THIS IS WHERE THE ACTION HAPPENS, WILL TAKE A LONG TIME
    // TO GET FIRST CONTOUR, THEN DO THE REST VERY FAST
    // ***********************************************************************
    TotalConts =  DoGetContour(streamFn, 10);



    printf("\n\n\tExtracted %d Graphs as Contours\n", TotalConts);
  
}




// *********************************************************************************
// *********************************************************************************
Int_t DoGetContour(TGraph2D *streamFn, Int_t nContours)
{  
    TList* contour      = NULL;  // List of Graphs for one contour     

    Int_t nTotalConts = 0;       // total # of graphs in all contours

    Int_t i             = 0;     // contour counter 
    Double_t level      = 0;     // current level
    Double_t max = 0, min = 0;   
    Double_t spacing = 0;

    if (streamFn->GetN() < 1)
      return -1 ; // Points haven't been properly calculated
    
    max = streamFn->GetZmax();
    min = streamFn->GetZmin();

    spacing = (max - min)/(nContours);

    // Get Contours from stream function
    for(i = 0; i < nContours; i++){       
       level = (max - 0.5*spacing) - i*spacing;

       printf("\n\n******************************\n   Level = %f    \n", level);
       printf("...GETTING CONTOUR LIST\n" );

       // Get list of Graphs making up this contour level
       contour = streamFn->GetContourList(level);

       printf("\tThis Level Has %d Contours\n", contour->GetSize(););
       printf("******************************\n");

       nTotalConts += contour->GetSize();
    }  // END FOR
       
    return nTotalConts;  // success
}



// *********************************************************************************
// *********************************************************************************
Double_t SawTooth(Double_t x, Double_t WaveLen){

/* This function is specific to a sawtooth function with period
  WaveLen symetric about x = 0 and with amplitude = 1. Each segment
  is 1/4 of the wavelength.

           |  
      /\   |
     /  \  |
    /    \ |
   /      \          
  /--------\--------/------------
           |\      /
           | \    / 
           |  \  /
  	   |   \/

*/
    Double_t y;
    if ( (x < -WaveLen/2) || (x > WaveLen/2))
      y = -99999999; // ERROR X OUT OF BOUNDS, Not elegant


    if (x <= -WaveLen/4){
	y = x + 2.0;
    } 
    else if ( (x > -WaveLen/4) && (x <= WaveLen/4)){
	y = -x ;
    } 
    else if( ( x > WaveLen/4) && (x <= WaveLen/2)){
	y = x - 2.0;
    }

    return y;
}