Object3d (xfig example)

Plotting in 3D

% -*- mode: slang; mode: fold -*-
require("isisscripts");
require("lbscripts");
()=evalfile("/home/falkner/work/lightbending/lbscripts.git/modules/ptrace.sl");

%%%%% GLOBAL PARAMETER: =======================================================

xfig_set_latex_preamble("\renewcommand{\vec}[1]{\ensuremath{\mathbf{#1}}}"R);

%%% CAMERA: modified Spherical angles
variable cam1 = 20*PI/180.; % cam1 = PI/2 -  thetaa
variable cam2 =  0*PI/180.; % cam2 = PI - phi


%%% LIGHT BENDING:
variable Mns = 1.44*2e30;
%variable Rns = 10e3; % Physical radius of the neutron star
variable rs = 2*Const_G*Mns/sqr(Const_c)*1e1; % Schwarzschild radius
variable r = 15e3; % radius of photon emission
variable u = rs/r;

%%% GEOMETRY:
variable R = 1; % radius of the sphere [BUG HERE FOR R /= 1]

variable i = 60*PI/180.; % Inclination of the observer

variable theta =  45*PI/180.; % Spot position
variable phi   = 220*PI/180.; % Spot position

%%% PLOT DIMENSION
variable size = 12;
variable xmin = -1.00*R;
variable xmax =  2.20*R;
variable ymin = -0.55*R;
variable ymax =  2.10*R;

%%% PLOT COLORS
variable col_cbg = ["#88CCEE",
		    "#999933",
		    "#AA4499",
		    "#332288"
		   ];
% variable photon1 = col_cbg[2];
% variable photon2 = col_cbg[3];
% variable obsplane = col_cbg[1];
variable photon1 = CB_COLOR_SCHEME[1];
variable photon2 = CB_COLOR_SCHEME[2];
variable obsplane = "#505050";


%%% PHYSIKAL VECTORS ==========================================================

%%% CAMERA: Basis vectors
variable c1 = vector_rotate( vector(1,0,0), vector(1,0,0), cam1 );
variable c2 = vector_rotate( vector(0,1,0), vector(1,0,0), cam1 );
variable c3 = vector_rotate( vector(0,0,1), vector(1,0,0), cam1 );
c1 = vector_rotate( c1, vector(0,0,1), cam2 );
c2 = vector_rotate( c2, vector(0,0,1), cam2 );
c3 = vector_rotate( c3, vector(0,0,1), cam2 );

%%% AXIS
variable ex = vector( R, 0, 0 );
variable ey = vector( 0, R, 0 );
variable ez = vector( 0, 0, R );

%%% OBSERVER
variable k = vector( sin(i), 0, cos(i) );
%%% SPOT
variable n  = vector( 1, phi, theta ; sph );

%%% LIGHT BENDING CALCULATIONS
%%% Attention: simple calculations based on beloborodov 2002, i.e., it
%%% does not work for all settings (where alpha > 90°)

variable lbpar = struct{
  nmr = 10,
  nma = 100,
  rmin = 0.95*r/rs,
  rmax = 1.05*r/rs,
  lbmeth = "belob",
  tdelay=0,
};
variable lbmap = lb_psir2b_map( lbpar );


%% Solution 1:
variable psi   = acos( dotprod(k,n) ); 
variable alpha = acos( 1 - ( 1-cos(psi) )*( 1-u ) );
alpha = lb_interp( lbmap, r/rs, psi, struct{ a, b } ).a;

variable b     = R*sin(alpha)/sqrt(1-u);
b = lb_interp( lbmap, r/rs, psi, struct{ a, b } ).b*rs/r*R;

variable crho  = ( ( k^(n^k) ) * ey )/abs( k^(n^k) );
variable srho  = abs( ( k^(n^k) ) ^ ey )/abs( k^(n^k) );

variable k0 = unit_vector(vector_rotate( n, crossprod(n,k), alpha ));

%% Solution 2:
%%   1 - cos(alpha) = ( 1 - cos(psi) )*(1-u)
%%   b = R*sin(alpha) / sqrt(1-u)
variable _psi   = 2*PI - psi;
variable _alpha = acos( 1 - ( 1-cos(_psi) )*( 1-u ) );;
_alpha = lb_interp( lbmap, r/rs, _psi, struct{ a, b } ).a;


variable fac = rs/r*R;
variable _b     = R*sin(_alpha)/sqrt(1-u);
_b = lb_interp( lbmap, r/rs, _psi, struct{ a, b } ).b*fac;
variable _rp = lb_p_b( _b/fac )*fac;
variable _pp = lb_psiper_u_belob_( fac/_rp );

variable _k0 = unit_vector(vector_rotate( n, crossprod(n,k), -_alpha ));




%%% FUNCTIONS =================================================================
%%% Must be defined here, as functions require some global variables!

% Projection of vectors onto the camera plane
define P( v ){ %{{{
  variable vec = @v;
  if( typeof(vec)== Array_Type ){
    vec = struct_array_2_struct_of_arrays( v );
  }
  variable pro = vector_change_basis( vec, c1, c2, c3 );
  return pro;
}
%}}}
% Return x,z component of Projection
define PP( v ){ %{{{
  variable p = P(v);
  return p.x, p.z;
}
%}}}
% GREAT CIRCLE between 2 Vector
define GC( a, b ){ %{{{
  
  variable N = qualifier("N",100);
  variable R = qualifier_exists("R")? 1:0;
  
  variable na = unit_vector(a);
  variable nb = unit_vector(b);

  variable rotangle = [0:acos(dotprod(na,nb))-R*2*PI:#N];
  variable rotvec   = unit_vector( crossprod( a, b ) );

  variable gc = vector_rotate( na, rotvec, rotangle );

  return gc;
}
%}}}
% PHOTON TRAJECTORY
define T(p){ %{{{
  variable B = b/R*r;
  variable X = rs*( 1 - cos(p))/(2*( 1 + cos(p)));
  variable trace = sqrt( sqr(X) + sqr(B/sin(p)) ) - X;
  return trace/r*R;
}
%}}}
define PT(){ %{{{
  variable brs = b/R*r/rs;
  variable pt = lb_ptrace_belob( brs ; n=N, pmax=psi );
  
  variable v = vector( pt.r*rs/r*R, [0:0:#N], -pt.p ; sph );
  v = vector_rotate( v, ey, i );
  v = vector_rotate( v, k, -atan2(crho,srho) );

  return v;  
}
%}}}
define _PT(){ %{{{
  variable brs = _b/R*r/rs;

  variable pt = lb_ptrace_belob( brs ; n=N, pmax=2*PI-psi );
  
  variable v = vector( pt.r*rs/r*R, [0:0:#N], -pt.p ; sph );
  v = vector_rotate( v, ey, i );
  v = vector_rotate( v, k, PI-atan2(crho,srho) );

  return v;  
}
%}}}

  



%%% PREPARE PLOT ==============================================================
variable s = 1.5;
variable ls= s*1.06;
variable N = 300;

% arrow head
variable ARROW = xfig_create_arrow(;arrow_type=2,
				   arrow_style=1,
				   arrow_width=8,
				   arrow_heigth=64,
				   arrow_thickness=2
				  );
variable AXIS = @ARROW;
AXIS.arrow_style=0;
AXIS.arrow_type=0;

% circle
variable xycirc = vector( [R:R:#N], [0:2*PI:#N], [PI/2.:PI/2.:#N]; sph );
variable xzcircn = vector( [R:R:#N], [0:0:#N], PI/2.*[-1:1:#N]; sph );
variable xzcircs = vector( [R:R:#N], [0:0:#N], PI/2.*[1:3:#N]; sph );

% n help lines
variable ncirc  = vector( [R:R:#N], [phi:phi:#N], [0:PI/2:#N]; sph );
variable nline  = vector( [0:R:#N], [phi:phi:#N], [PI/2:PI/2:#N];sph);

% kn circ
variable kncirc = R*GC( unit_vector(n), unit_vector(k) );
variable nkcirc = R*GC( unit_vector(n), unit_vector(k) ; R );

%%% PLOT ======================================================================
variable xp = xfig_plot_new(size,(ymax-ymin)/(xmax-xmin)*size);
xp.world( xmin, xmax, ymin, ymax ; padx=0.0, pady=0.0 );
xp.axis(;off);

% axis
xp.plot( PP([0,s]*R*ex); forward_arrow=AXIS );
xp.xylabel( PP(ls*R*ex), "$x$"R  );

xp.plot( PP([0,s]*ey) );
%xp.xylabel( PP(ls*ey), "$y$"R  );

xp.plot( PP([0,s]*R*ez); forward_arrow=AXIS );
xp.xylabel( PP(ls*ez), "$z$"R  );

% circles
xp.plot( PP(xzcircn) ; line=0, width=2 );
xp.plot( PP(xzcircs) ; line=0, width=2 );
xp.plot( PP(xycirc)  ; line=0, width=2 );

% rotation
variable rot = 0.8*s*R*ez+vector( 0.15*[R:R:#N], PI/180*[-60:265:#N], [PI/2.:PI/2.:#N]; sph );
xp.plot( PP(rot); forward_arrow=AXIS );

% observer
xp.plot( PP([0,1]*R*k) ; line=1 );
xp.plot( PP([1,s]*R*k) ; line=0 );
xp.xylabel( PP(ls*k), "$\vec{k}$"R, 1.5, 1.5 );
xp.plot( PP([1,s]*R*k); forward_arrow=ARROW );

% observer plane
variable deltap = 1.8*R;
variable oplane = 1.2*R*vector( [-1, 1, 1,-1,-1],
			  [-1,-1, 1, 1,-1],
			  [ 0, 0, 0, 0, 0]
			);
oplane.z += deltap;
oplane = vector_rotate( oplane, ey, i );
xp.plot( PP(oplane); depth=40, width=3, color=obsplane );
%variable oplab = 1.8*R*vector( -.01, -.95, 0 );
variable oplab = 1.8*R*vector( .60, 0, 0 );
oplab.z += deltap;
oplab = vector_rotate( oplab, ey, i );
%xp.xylabel( PP(oplab), "observer plane"R ; rotate=-0.97*i*180/PI, color=obsplane  );
xp.xylabel( PP(oplab), "{\footnotesize observer plane}"R ; rotate=90, color=obsplane  );
variable ocirc = vector( [R:R:#N],[0:2*PI:#N],PI/2.*[1:1:#N]; sph );
ocirc.z += deltap;
ocirc = vector_rotate( ocirc, ey, i );
xp.plot( PP(ocirc); color=obsplane );
variable A = 1.15*s*R*vector( [0,0],[0,1],[0,0] );
A.z += deltap;
A = vector_rotate( A, ey, i );
xp.plot( PP(A) ; forward_arrow=AXIS, color=obsplane );
xp.xylabel( PP(struct_of_arrays_2_struct_array(A)[-1]), "$A$"R, -.5, 1.2;color=obsplane);
variable B = s*R*vector( [0,-1],[0,0],[0,0] );
B.z += deltap;
B = vector_rotate( B, ey, i );
xp.plot( PP(B) ; forward_arrow=AXIS,color=obsplane );
xp.xylabel( PP(struct_of_arrays_2_struct_array(B)[-1]), "$B$"R, -.4, 2.2;color=obsplane );


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% IMPACT A & B
%%%%%%%%%%%%%%%%

% Line in OP from center to impact point 1
variable IB = b*vector( [0,1], PI/2+[0,PI/2-atan2(crho,srho)], PI/2*[1,1] ; sph );
variable IB2 = @IB;
IB.z += deltap;
IB = vector_rotate( IB, ey, i );
xp.plot( PP(IB) ; line=1,color=obsplane );
% Projected circle at the impact point 1
variable IBcirc = vector( 0.05*[R:R:#N],[0:2*PI:#N],PI/2.*[1:1:#N]; sph );
IBcirc.z += deltap;
IBcirc = vector_rotate( IBcirc, ey, i );
IBcirc += struct_of_arrays_2_struct_array(IB)[-1] - struct_of_arrays_2_struct_array(IB)[0];
IBcirc = struct_of_arrays_2_struct_array(IBcirc);
%xp.shade_region( PP(IBcirc) ; fillcolor=photon1,width=0);
xp.plot( PP(struct_of_arrays_2_struct_array(IB)[-1]) ; sym="circle", color=photon1,fill, size=0.3*R);

% Arrow between k and photon1 tratecory indicating b
IB2.z += 0.7*deltap;
IB2 = vector_rotate( IB2, ey, i );
xp.plot( PP(IB2) ; line=0 , forward_arrow, backward_arrow, depth=150 );
xp.xylabel( PP(vector( mean(IB2.x), mean(IB2.y), mean(IB2.z) )), "$b$"R, -.4, -.4 );

% Line in OP from center to impact point 2
variable _IB = -_b*vector( [0,1], PI/2+[0,PI/2-atan2(crho,srho)], PI/2*[1,1] ; sph );
_IB.z += deltap;
_IB = vector_rotate( _IB, ey, i );
xp.plot( PP(_IB) ; line=1,color=obsplane );
% Projected circle at the impact point 1
variable _IBcirc = vector( 0.05*[R:R:#N],[0:2*PI:#N],PI/2.*[1:1:#N]; sph );
_IBcirc.z += deltap;
_IBcirc = vector_rotate( _IBcirc, ey, i );
_IBcirc += struct_of_arrays_2_struct_array(_IB)[-1] - struct_of_arrays_2_struct_array(_IB)[0];
_IBcirc = struct_of_arrays_2_struct_array(_IBcirc);
%xp.shade_region( PP(_IBcirc) ; fillcolor=photon2,width=0);
xp.plot( PP(struct_of_arrays_2_struct_array(_IB)[-1]) ; sym="circle", color=photon2,fill, size=0.3*R);


% RHO
variable rhoc = 0.3*R*vector( [1:1:#N], PI/2+[0:PI/2-atan2(crho,srho):#N], PI/2*[1:1:#N] ; sph );
rhoc.z += deltap;
rhoc = vector_rotate( rhoc, ey, i );
xp.plot( PP(rhoc) ; line=0,color=obsplane );
xp.xylabel( PP(vector( mean(rhoc.x), mean(rhoc.y), mean(rhoc.z) )), "$\rho$"R, -.4, -.8;color=obsplane );


xp.plot( PP([R,deltap]*k) ; line=2, width=1 );



% inclination
variable ic = GC( ez, k );
xp.plot( PP(0.18*ic) );
xp.xylabel( PP(0.12*vector( mean(ic.x), mean(ic.y), mean(ic.z) )), "$i$"R );

% PSI
variable pc = GC( n, k );
xp.plot( PP(0.42*pc) ; color=photon1, depth=50);
xp.xylabel( PP(0.35*vector( mean(pc.x), mean(pc.y), mean(pc.z) )), "$\varPsi$"R, -.5, 0; color=photon1 );

variable _pc = GC( n, k ; R);
xp.plot( PP(0.33*_pc) ; color=photon2, depth=150 );
xp.xylabel( PP(1.*vector( mean(_pc.x), mean(_pc.y), mean(_pc.z) )), "$\varPsi^*$"R, -1.8, 0 ; color=photon2);

% Periastron
variable _per = _rp*vector_rotate( unit_vector(n), crossprod(k,n), _psi-_pp ); 
variable _pper = GC( _per, k);
xp.plot( PP([0,1]*_per) ; color=photon2 );
xp.xylabel( PP(0.6*_per), "$r_\mathrm{p}$"R, .7, .0 ; color=photon2);
xp.plot( PP(0.29*_pper) ; color=photon2, depth=150 );
xp.xylabel( PP(0.25*vector( mean(_pper.x), mean(_pper.y), mean(_pper.z) )), "$\varPsi_\mathrm{p}$"R, .25, .0 ; color=photon2);


%%% Emission point: R, n & help lines
xp.plot( PP([0,1]*R*n) ; line=1 );
xp.xylabel( PP(0.5*R*n), "$R$"R, .5, .5 );
xp.plot( PP([1,1.2*s]*R*n) );
xp.xylabel( PP(1.2*ls*R*n), "$\vec{n}$"R, -1,0 );
xp.plot( PP([1,1.2*s]*R*n); forward_arrow=ARROW );

xp.plot( PP(ncirc) ; line=2 );
xp.plot( PP(nline) ; line=2 );
xp.plot( PP(kncirc); line=1 );
xp.plot( PP(nkcirc); line=1 );

variable Scirc = vector( 0.05*[R:R:#N],[0:2*PI:#N],PI/2.*[1:1:#N]; sph );
Scirc.z += R;
Scirc = vector_rotate( Scirc, unit_vector(crossprod(ez,n)), theta );
Scirc = struct_of_arrays_2_struct_array(Scirc);
%xp.shade_region( PP(Scirc) ; fillcolor="black",width=0);
xp.plot( PP(R*n)  ; sym="circle", fill, size=0.3*R, depth=40 );

% phi
variable phic = vector( 0.13*[R:R:#N], [0:phi:#N], [PI/2.:PI/2:#N]; sph );
xp.plot( PP(phic) );
xp.xylabel( PP(vector( mean(phic.x), mean(phic.y), mean(phic.z) )), "$\phi$"R, 1.2, -.1 );

% theta
variable tc = vector( 0.6*[R:R:#N], [phi:phi:#N], [0:theta:#N]; sph );
xp.plot( PP(tc) );
xp.xylabel( PP(0.8*vector( mean(tc.x), mean(tc.y), mean(tc.z) )), "$\theta$"R, -.1,-.5 );


% photon
xp.plot( PP(R*[n,n+(s-1)*k0]) ; line=0 , depth=40);
xp.xylabel( PP(R*(n+(ls-1)*k0)), "$\vec{k}_0$"R, 1.6, .6; depth=40 );
xp.plot( PP(R*[n,n+(s-1)*k0]); forward_arrow=ARROW, depth=40);

% alpha
variable ac = GC( k0, n );
xp.plot( PP(0.2*ac+R*n) );
xp.xylabel( PP(R*n+0.13*vector( mean(ac.x), mean(ac.y), mean(ac.z) )), "$\alpha$"R );

% trace
variable t1 = PT();
variable t1ind = where( vector_norm(t1) <= abs(vector(IB.x[-1],IB.y[-1],IB.z[-1])) );
xp.plot( PP(struct_filter( t1, t1ind ; copy ) ) ; color=photon1, depth=50, width=3 );
xp.plot( PP(struct_filter( t1, complement([0:length(t1ind)-1],t1ind) ; copy )); color=photon1, depth=30, line=1 );
variable t2 = _PT();
variable t2ind = where( vector_norm(t2) <= abs(vector(_IB.x[-1],_IB.y[-1],_IB.z[-1])) );
xp.plot( PP(struct_filter( t2, t2ind ; copy ) ) ; color=photon2, line=1,depth=250, width=3 );
xp.plot( PP(struct_filter( t2, complement([0:length(t2ind)-1],t2ind) ; copy ) ); color=photon2, depth=30, line=1 );


%%% FILL OBSERVER PLANE
variable fill1 = vector([kncirc.x,IB.x[[0,1]],struct_filter( t1, t1ind ; copy ).x],
			[kncirc.y,IB.y[[0,1]],struct_filter( t1, t1ind ; copy ).y],
			[kncirc.z,IB.z[[0,1]],struct_filter( t1, t1ind ; copy ).z]
		       );
variable _fill1 = vector([0,IB.x[[0,1]],struct_filter( t1, t1ind ; copy ).x],
			 [0,IB.y[[0,1]],struct_filter( t1, t1ind ; copy ).y],
			 [0,IB.z[[0,1]],struct_filter( t1, t1ind ; copy ).z]
			);
variable fill2 = vector([nkcirc.x,_IB.x[[0,1]],struct_filter( t2, t2ind ; copy ).x],
			[nkcirc.y,_IB.y[[0,1]],struct_filter( t2, t2ind ; copy ).y],
			[nkcirc.z,_IB.z[[0,1]],struct_filter( t2, t2ind ; copy ).z]
		       );
variable _fill2 = vector([0,_IB.x[[0,1]],struct_filter( t2, t2ind ; copy ).x],
			 [0,_IB.y[[0,1]],struct_filter( t2, t2ind ; copy ).y],
			 [0,_IB.z[[0,1]],struct_filter( t2, t2ind ; copy ).z]
			);
xp.shade_region( PP(_fill1); fillcolor=photon1, fill=36,depth=250, width=0 );
xp.shade_region( PP(_fill2); fillcolor=photon2, fill=36,depth=250, width=0 );
		  

xp.render("geo.pdf");
xp.render("geo.png");