share/tensor/kruskal.dem

   1 /* Copyright (C) 2008 Viktor T. Toth <http://www.vttoth.com/>
   2  *
   3  * This program is free software; you can redistribute it and/or
   4  * modify it under the terms of the GNU General Public License as
   5  * published by the Free Software Foundation; either version 2 of
   6  * the License, or (at your option) any later version.
   7  *
   8  * This program is distributed in the hope that it will be
   9  * useful, but WITHOUT ANY WARRANTY; without even the implied
  10  * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11  * PURPOSE.  See the GNU General Public License for more details.
  12  *
  13  * The Kruskal extension of the Schwarzschild metric.
  14  *
  15  */
  16 kill(all);
  17
  18 ("Computing Kruskal's maximal extension of the Schwarzschild metric" )$
  19
  20 ("We begin with redefinining the derivative of Lambert's W function.")$
  21 gradef(lambert_w(x),lambert_w(x)/(1+lambert_w(x))/x);
  22 ("We want ctensor to perform rational factorization.")$
  23 ratfac:true;
  24 load(ctensor);
  25 ct_coords:[T,R,u,v];
  26 ("This is Kruskal's metric in explicit form")$
  27 r:2*(lambert_w((R-T)*(R+T)*exp(-1))+1)*G*M;
  28 lg:ident(4);
  29 lg[1,1]:-lg[2,2]:-16*G^2*M^2/(-T^2+R^2)*lambert_w(exp(-1)*(-T^2+R^2))/(lambert_w(exp(-1)*(-T^2+R^2))+1);
  30 lg[3,3]:-r^2;
  31 lg[4,4]:-r^2*sin(u)^2;
  32 ric:zeromatrix(4,4);
  33 cmetric();
  34 christof(mcs);
  35 ricci(false);
  36 ("Let us verify that this is a vacuum metric (Ricci tensor is zero)")$
  37 radcan(ric);
  38 ("We now wish to express the metric coefficients as functions of r")$
  39 kill(r);
  40 ("This is t: we can now plot the tt component of the metric")$
  41 block([title: "The Kruskal-Szekeres metric coefficient for GM=1"],
  42 /*plot3d([r*cos(u),r*sin(u),ev(subst((1-r/(2*G*M))*exp(r/(2*G*M))+R^2,T^2,ug[1,1]),G=1,M=1)],[r,0.01,5],[u,-%pi,%pi],['grid,50,15]));*/
  43 /*plot3d([r*cos(u),r*sin(u),ev(subst((1-r/(2*G*M))*exp(r/(2*G*M))+R^2,T^2,ug[1,1]),G=1,M=1)],[r,0.01,5],[u,-%pi,%pi],['grid,7,7],[plot_format,gnuplot],[gnuplot_term,dumb],[gnuplot_dumb_term_command, "set term dumb 132 50"]));*/
  44 plot2d(ev(subst((1-r/(2*G*M))*exp(r/(2*G*M))+R^2,T^2,ug[1,1]),G=1,M=1),[r,0,5],[plot_format,gnuplot],[gnuplot_term,dumb],[yx_ratio,0.25]));
  45 /* End of demo -- comment line needed by MAXIMA to resume demo menu */