//###########################################################################
//
// Usage: 
//  - launch ROOT (you need version 3.05.06)
//  - at root prompt, do:
//
//      root[0] gSystem->Load("libGeom");
//      root[1] .x polar.C
//
//  - you can rotate the view with the mouse and at any time repeat the
//    animation :
//      root[3] Animate(nframes); // default is 200 frames
//    or just draw the fields:
//      root[4] DrawEH(Bool_t e, Bool_t h) // by default both true - see both
//
//###########################################################################


Double_t lambda = 1; // should be ~1E-9, but geometry cannot handle this value
Double_t c = 1E8;    // wave speed in vacuum
Double_t ehAspectRatio = 0.5; // put 1 if you want to have the same amplitude
Double_t tmin = 0;
Double_t tmax = 10*lambda/c;  // time range
Int_t ie,ih;
Double_t fx = 1.0;  // bigger fx = bigger sines

void polar()
{
// create a simple geometry with a quarter wave plate disposed at z=5*lambda,
// and having a thickness of lambda/4
   TGeoManager *geom = new TGeoManager("NA59", "quarter wave plate");
   // dummy materials and media - no tracking physics
   TGeoMaterial *mat = new TGeoMaterial("Vacuum",0,0,0);
   TGeoMedium   *med = new TGeoMedium("Vacuum",1,mat);
   TGeoVolume *top = geom->MakeBox("TOP", med, 30.*lambda, 30.*lambda, 30.*lambda);
   top->SetVisibility(kFALSE);
   geom->SetTopVolume(top);
   TGeoVolume *quarter = geom->MakeBox("QUARTER", med, 2*lambda, 2*lambda, 0.25*lambda);
   quarter->SetLineColor(4);
   quarter->SetLineWidth(2);
   // position the quarter at 5*lambda on Z axis
   top->AddNode(quarter, 1, new TGeoTranslation(0,0,5*lambda));
   geom->CloseGeometry();
   
   // now create the electric and magnetic field "tracks"
   ie = geom->AddTrack(1,1);
   ih = geom->AddTrack(1,2);
   geom->SetPdgName(1,"E");
   geom->SetPdgName(2,"H");
   TVirtualGeoTrack *trackE = geom->GetTrack(ie);
   trackE->SetLineColor(kRed);
   trackE->SetMarkerColor(kRed);
   TVirtualGeoTrack *trackH = geom->GetTrack(ih);
   trackH->SetLineColor(kBlue);
   trackH->SetMarkerColor(kBlue);
   Int_t istep;
   // compute 200 points (Ex, Ey) -> (x,y,z) according to:
   //   z = c*t
   //   kz-wt = (kc-w)*t = (kc-w)*z/c = 2*pi*z/lambda
   //   Ex(z) = E0*cos(2*pi*z/lambda)
   //   Ey(z) = E0*sin(2*pi*z/lambda)
   //
   // For visualization reasons, we need several real points (x,y,z) at different
   // moments t, expressed in the geometry coordinate system. We therefore will
   // make an arbitrary conversion from the amplitude of the wave at an amplitude
   // on x or y axis, equal to 0.5*lambda:
   //   x(z) = fx*lambda*cos(2*pi*z/lambda);
   //   y(z) = fy*lambda*sin(2*pi*z/lambda);
   //   t(z) = z/c
   // We then sample five wave lengths before and after the quarter wave layer in
   // 100 points each. The full time range is (0,10.*lambda/c)
   Double_t z = 0;
   Double_t dz = lambda/20.;
   Double_t point[4]; // (x,y,z,t)
   Double_t dummy;
   for (istep=0; istep<200; istep++) {
      z = istep*dz;
      point[2] = z;   // z position
      point[3] = z/c; // time
      if (z<5*lambda) point[1]=0;  // linear polarization on X before quarter
      else point[1] = fx*lambda*TMath::Cos(2.*TMath::Pi()*z/lambda); // y amplitude
      point[0] = fx*lambda*TMath::Sin(2.*TMath::Pi()*z/lambda); // x amplitude
      trackE->AddPoint(point[0],point[1],point[2],point[3]);
      // invert x-y for H
      dummy = point[0];
      point[0] = point[1];
      point[1] = dummy;
      point[0] *= ehAspectRatio;
      point[1] *= ehAspectRatio;
      trackH->AddPoint(point[0],point[1],point[2],point[3]);
   } 
   // now draw the geometry  
   top->Draw();
//   geom->DefaultAngles();
   TView *view = gPad->GetView();
   view->Top();
   view->ShowAxis();
   Animate(50, "/*/S");
   DrawEH();
   new TBrowser();
}

void DrawEH(Bool_t e=kTRUE, Bool_t h=kTRUE)
{
   gGeoManager->GetTopVolume()->Draw();
   gPad->GetView()->ShowAxis();
   gGeoManager->SetTminTmax();
   if (e && h) {
      gGeoManager->DrawTracks();
      return;
   }
   if (e) {
      TVirtualGeoTrack *trackE = gGeoManager->GetTrack(ie);
      trackE->Draw();
      return;
   }   
   if (h) {
      TVirtualGeoTrack *trackH = gGeoManager->GetTrack(ih);
      trackH->Draw();
   }   
}

void Animate(Int_t nframes=2000, Option_t *option="/*")
{
      gGeoManager->AnimateTracks(tmin,tmax, nframes, option);
      gGeoManager->SetTminTmax();
}