Lightbending ptrace (xfig example)

Relativistic photon trajectory

require("/home/falkner/Public/lightbending.git/lbscripts");

define lb_plot_ptracesketch( R, b ){

  variable size = qualifier("size",16);
  variable rmax = qualifier("rmax",10.);
  variable index = qualifier("index","");
  if( index !="" )
    index = string( "\\ensuremath{_"+index+"}" );

  variable dx = qualifier("dx",1.);
  
  variable a1 = lb_a_ub( 1./R, b );
  variable a2 = PI-a1;
  variable p1 = lb_psi_ua_belob( 1./R, a1 );
  variable p2 = lb_psi_ua_belob( 1./R, a2 ; b=b );
  variable p = lb_p_b( b );
  variable psip = lb_psi_ua_belob( 1./p, PI/2. ; b=b );

  variable pon = 0;
  if( qualifier_exists("periastron") )
    pon = 1;
  variable pmax = p1;
  if(pon)
    pmax = p2;
  
  variable trace = lb_ptrace_belob( b; pmax=pmax, rmax=rmax);

  variable xmax = abs(max(trace.r)*cos(min(trace.p)));
  variable ymax = abs(max(trace.r)*sin(min(trace.p)));
  variable tmp,dr;

  
  variable xf;
  if( qualifier_exists("addto") )
    xf=qualifier("addto");
  else{
    xf = xfig_plot_new(size,size);
    xf.world(-1.2*xmax,1.2*xmax,-1.2*xmax,1.2*xmax;xpad,ypad);
    xf.axis(;off);
  }
  
  %Circle at R
  tmp = [0:2*PI:#100];
  xf.plot( R*cos(tmp), R*sin(tmp) ; line=2 );
  tmp = 2.5/R*7./4.*PI;
  xf.plot( [0,R*cos(tmp)], [0,R*sin(tmp)] ; line=2 );
  xf.xylabel(.5*R*cos(0.93*tmp), .5*R*sin(0.93*tmp) , "$R$"+index );
  
  %trace
  xf.plot( trace.r*cos(trace.p), trace.r*sin(trace.p) ; width=4, color="red", depth=20);
  xf.xylabel(dx*xmax,b,"photon"+index,.5,-.5);
  
  % impact parameter
  xf.plot( dx*[0.8*xmax, 0.8*xmax],[0,b] );
  xf.xylabel( dx*0.8*xmax, b/2., "$b$"+index, .5,0);

  ifnot( qualifier_exists("addto") ){
    % observer sky
    xf.plot( [1.1*xmax,1.1*xmax],[-1.2*ymax,1.2*ymax] ; line=3 );
    xf.xylabel( 1.12*xmax, 0, "{\\tiny observer sky at infinity}",-.5,0 ; rotate=90 );
    
    % line of sight
    xf.plot( [0,xmax],[0,0] );
    xf.xylabel( .5*(xmax+R),-.05*R,"\\tiny line of sight",.0,.5);
  }
  ifnot( qualifier_exists("nopsi") ){
    % Psi R
    xf.plot( [0,1.5*R*cos(p1)],[0,1.5*R*sin(p1)] ; depth=19 );
    tmp = [0:p1:#100];
    dr = .45;
    xf.plot( dr*R*cos(tmp),dr*R*sin(tmp) );
    dr += .2;
    xf.xylabel( dr*R*cos(p1/2.), dr*R*sin(p1/2.), "$\\Psi_{R"+index+"}$");
  }
  
  % alpha
  if( qualifier_exists("alpha") ){
    tmp = [p1-a1:p1:#100];
    dr=.4;
    xf.plot( R*cos(p1)+dr*R*cos(tmp), R*sin(p1)+dr*R*sin(tmp) ;depth=19);
    xf.xylabel( R*cos(p1)+.7*dr*R*cos(p1-a1/2.), R*sin(p1)+.7*dr*R*sin(p1-a1/2.),
		"$\\alpha$"+index );
    dr = 1.;
    xf.plot( R*cos(p1)+[0,dr*R*cos(p1-a1)],R*sin(p1)+[0,dr*R*sin(p1-a1)] ; line=1,depth=19);
  }

  if(pon){
    ifnot( qualifier_exists("nopsi") ){
      % Psi R prime
      xf.plot( [0,1.5*R*cos(p2)],[0,1.5*R*sin(p2)] ; depth=19);
      tmp = [0:p2:#100];
      dr = .15;
      xf.plot( dr*R*cos(tmp),dr*R*sin(tmp) );
      dr += .3;
      xf.xylabel( dr*R*cos(.5*(psip+p2)), dr*R*sin(.5*(psip+p2)), "$\\Psi^*_{R"+index+"}$");
      
      % periastron
      xf.plot( [0,p*cos(psip)],[0,p*sin(psip)] );
      tmp = [0:psip:#100];
      dr = .3;
      xf.plot( dr*R*cos(tmp), dr*R*sin(tmp) );
      dr += .2;
      xf.xylabel( dr*R*cos(.5*(psip+p1)), dr*R*sin(.5*(psip+p1)), "$\\Psi_{p"+index+"}$");
    }
    if( qualifier_exists("alpha") ){
      % alpha prime
      tmp = [p2-a2:p2:#100];
      dr = .4;
      xf.xylabel( R*cos(p2)+.6*dr*R*cos(p2-a2/2.), R*sin(p2)+.6*dr*R*sin(p2-a2/2.),
		  "$\\alpha^*$"+index );
      xf.plot( R*cos(p2)+dr*R*cos(tmp), R*sin(p2)+dr*R*sin(tmp) ;depth=19);
      dr = 1.5;
      xf.plot( R*cos(p2)+[0,dr*R*cos(p2-a2)],R*sin(p2)+[0,dr*R*sin(p2-a2)] ; line=1,depth=19);

      if(qualifier_exists("symm")){
	% alpha prime at alpha
	tmp = [p1:p1+a2:#100];
	dr = .4;
	xf.xylabel( R*cos(p1)+.6*dr*R*cos(p1+a2/2.), R*sin(p1)+.6*dr*R*sin(p1+a2/2.),
			  "$\\alpha^*$"+index );
	xf.plot( R*cos(p1)+dr*R*cos(tmp), R*sin(p1)+dr*R*sin(tmp);depth=19 );
	dr = 1.;
	xf.plot( R*cos(p1)+[0,dr*R*cos(p1+a2)],R*sin(p1)+[0,dr*R*sin(p1+a2)] ; line=1,depth=19);
      }
    }
  }
  return xf;
}

variable R = 2.5;
variable b = 2.9;

variable xf;
xf = lb_plot_ptracesketch( R, b );
xf.render("ptrace1.pdf");

xf = lb_plot_ptracesketch( R, b ; periastron );
xf.render("ptrace2.pdf");

xf = lb_plot_ptracesketch( R, b ; alpha,periastron,symm );
xf.render("ptrace3.pdf");

variable R2 = 1.5*R;
variable b2 = 1.3*b;
variable p1 = lb_psi_ua_belob( 1./R,   lb_a_ub( 1./R, b )  );
variable p2 = lb_psi_ua_belob( 1./R2, lb_a_ub( 1./R2, b2 ) );
variable p3 = lb_psi_ua_belob( 1./R2, PI-lb_a_ub( 1./R2, b2 ));
xf = lb_plot_ptracesketch( R2, b2 ; nopsi,index="2",periastron,dx=0.9);
xf = lb_plot_ptracesketch( R, b ; nopsi, addto=@xf, index="1" );
xf.plot( R,0; sym=1, size=1,fill );
xf.plot( R*cos(p1),R*sin(p1); sym=1, size=1,fill,fillcolor="red",depth=18 );
xf.plot( R2,0; sym=0, size=4,fill );
xf.plot( R2*cos(p2),R2*sin(p2); sym=0, size=4,fill,fillcolor="red",depth=18 );
xf.plot( R2*cos(p3),R2*sin(p3); sym=0, size=4,fill,fillcolor="white",depth=18 );
xf.render("ptrace4.pdf");

exit;