void earth()
{
   auto color_rgb = [](double r, double g, double b)
   {
      auto col_index = TColor::GetFreeColorIndex();
      new TColor(col_index, r/256, g/256, b/256);
      return col_index;
   };
   auto const earthcol = color_rgb(79, 76, 176);
   auto const linecol = kGray;
   auto const polecol = kWhite;
   constexpr double bbox_size = 100; // bounding box size
   constexpr double R = 90; // earth radius < bbox_size
   constexpr double w = 0.15; // meridian/parallel width

   auto geom = new TGeoManager("earth", "");
   TGeoMedium *med = new TGeoMedium("earth_med", 1, new TGeoMaterial("earth_average", 26, 55, 5.513));
   TGeoVolume *top = geom->MakeBox("universe",med, bbox_size, bbox_size, bbox_size);
   geom->SetTopVolume(top);

   TGeoVolume *earth_vol = geom->MakeSphere("earth",med, 0, R);
   earth_vol->SetLineColor(earthcol);
   top->AddNode(earth_vol,1);

   // meridians
   TGeoVolume *meridian = geom->MakeTube("meridian",med, R, R + w, w);
   meridian->SetLineColor(linecol);
   for (auto i = 0; i < 18; ++i) {
      double phi = 10 * i; // degrees
      auto rot = new TGeoRotation();
      rot->RotateX(90);
      rot->RotateZ(phi);
      top->AddNode(meridian, i, rot);
   }
   
   // poles
   TGeoVolume *pole = geom->MakeTube("pole",med, 0, 4 * w, w);
   pole->SetLineColor(polecol);
   top->AddNode(pole, 0, new TGeoTranslation(0, 0, R + w));
   top->AddNode(pole, 1, new TGeoTranslation(0, 0, -R - w));

   // equator
   TGeoVolume *equator = geom->MakeTube("equator",med, R, R + w, 2 * w);
   equator->SetLineColor(polecol);
   top->AddNode(equator, 0);

   auto make_parallels = [&](double theta_deg, double wi, int col)
   {
     // make parallels at angle theta and 180 - theta
      double theta = theta_deg * TMath::DegToRad();
      double r = R * std::sin(theta);
      double z = R * std::cos(theta);
      double dz = wi * r / R;
      double dr = wi * z / R;
      double rmin1 = r + dr;
      double rmax1 = r + dr + wi;
      double rmin2 = r - dr;
      double rmax2 = r - dr + wi;
      TGeoVolume *parallel_t = geom->MakeCone("parallel",med, dz, rmin1, rmax1, rmin2, rmax2);
      parallel_t->SetLineColor(col);
      top->AddNode(parallel_t, 0, new TGeoTranslation(0, 0, z));
      TGeoVolume *parallel_b = geom->MakeCone("parallel",med, dz, rmin2, rmax2, rmin1, rmax1);
      parallel_b->SetLineColor(col);
      top->AddNode(parallel_b, 0, new TGeoTranslation(0, 0, -z));
   };

   for (auto i = 1; i < 9; ++i) {
      make_parallels(10 * i, w, linecol);
   }

   // polar circles
   make_parallels(33.3, 2 * w, polecol);

   geom->CloseGeometry();
   geom->SetNsegments(256);
   top->Draw("ogl");
   //top->Raytrace();
}