// TEST SCRIPT FOR: Getting Contours From TH2D
// Author: Josh de Bever
//         CSI Medical Physics Group
//         The University of Western Ontario
//         London, Ontario, Canada
//   Date: Oct. 22, 2004


OneContNoPlot(){
    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, 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);

   
    // 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]); 
    }

    // Create Histogram
    TH2D *HistStreamFn = new TH2D("HstreamFn", 
				  "Histogram Stream Function", 
				  nZsamples, z[0], z[nZsamples-1],
				  nPhiSamples, phi[0], phi[nPhiSamples-1]);
     
    printf("Input Functions Calculated\n");
     
    // Load Histogram Data
    for (Int_t i = 0; i < nZsamples; i++) {
    	for(Int_t j = 0; j < nPhiSamples; j++){
    	    HistStreamFn->SetBinContent(i,j, HofZ[i]*FofPhi[j]);
    	} 
    }

    gStyle->SetPalette(1); 
    c->Divide(2,1);

    c->cd(1);
    // Draw contours as Lines
    HistStreamFn->Draw("CONT1");

    c->cd(2);

    Double_t contours[4];
  
    contours[0] = -0.4;
     
    HistStreamFn->SetContour(1, contours); 
  
    // Draw contours as filled regions, and Save points
    HistStreamFn->Draw("CONT Z LIST");
    c->Update();    

    // Get Contours
    TObjArray *conts = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
    TList* contLevel = NULL;
    TList* contLevel = NULL;
    TGraph* curv = NULL;    
    Int_t nGraphs = 0;
    
   
    if (conts == NULL){
    	printf("***No Contours Were Extracted!\n");
    	TotalConts = 0;
    }     
    else{
    	TotalConts = conts->GetSize();
    }
    
    for(i = 0; i < TotalConts; i++){       
    	contLevel = (TList*)conts->At(i);
	printf("\nContour %d has %d Graphs\n", i, contLevel->GetSize());


    	printf("ContLevel[%d]:  Z = %f\n", i, contours[i]);

	// Get first graph from list on curves on this level
	curv = (TGraph*)contLevel->First();
	for(j = 0; j < contLevel->GetSize(); j++){
	    nGraphs ++;
	    
	    printf("\tGraph: %d  -- %d Elements\n", nGraphs,  curv->GetN());
	   	    
	    // Add graph to list of graphs
	    curv = (TGraph*)contLevel->After(curv);
	    
	}

    }

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

// *********************************************************************************
// *********************************************************************************
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 but no
                     // other way to return error

    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;
}