;+
; NAME:
; SRVADD
;
; AUTHOR:
; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
; craigm@lheamail.gsfc.nasa.gov
; UPDATED VERSIONs can be found on my WEB PAGE:
; http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
; Add velocity 3-vectors according to special relativity
;
; MAJOR TOPICS:
; Physics, Geometry
;
; CALLING SEQUENCE:
; U0 = SRVADD(U1, V)
;
; DESCRIPTION:
;
; The function SRVADD performs addition of velocity 3-vectors
; according to special relativity.
;
; Consider two inertial coordinate frames. Frame "0" is a "lab" or
; rest frame. Frame "1" is a "rocket" or moving frame, moving at
; velocity V with respect to the lab frame. The velocity V is
; allowed to be an arbitrary 3-vector.
;
; * An observer in the rocket frame sees a body moving at velocity U1.
;
; * An observer in the lab frame sees the same body moving at
; velocity U0.
;
; * This function solves for U0 in terms of V and U1.
;
; U1 and V are allowed to be 3xN arrays, which means more than one
; vector can be computed in a single call. If the dimensions of
; either U1 or V are 3x1, then it will be expanded to match the
; dimensions of the other vector. This simulates addition by a
; "scalar" vector. Because V can be a 3xN array, this means that
; multiple "rocket" frames can be computed at one time.
;
; NOTE: Velocities passed to SRVADD are measured as a *fraction of
; the speed of light*. Therefore, if the velocities are
; measured in some physical units, and CLIGHT is the speed of
; light in those same units, then the following statement:
;
; U0 = SRVADD(U1/CLIGHT, V/CLIGHT)*CLIGHT
;
; will compute the velocity U0, also in the same units.
;
;
; The formula for computing the velocity in the lab frame is:
;
; ( (1-1/GAMMA)*(U1 . VUNIT)*VUNIT + U1/GAMMA + V )
; U0 = -------------------------------------------------
; (1 - U1 . V)
;
; where
; GAMMA is the Lorentz factor = 1/SQRT(1 - |V|^2)
; VUNIT is the unit vector in the direction of V, = V/|V|
; "." is the vector dot product
;
; [ IDL notation is not strictly adhered to in this formula, for
; clarity of presentation. ]
;
;
; INPUTS:
;
; U1 - 3-vector or 3xN array, the velocity of a body as seen in the
; rocket frame (frame 1). The velocity is normalized such that
; the speed of light is 1.
;
;
; V - 3-vector or 3xN array, the velocity of the rocket frame as
; seen by an observer in the lab. The velocity is normalized
; such that the speed of light is 1.
;
; RETURNS:
;
; A 3xN array, containing the velocity of the body as seen in the
; lab frame. The velocity is normalized such that the speed of
; light is 1.
;
; KEYWORD PARAMETERS:
;
; CLASSICAL - if set, then classical velocity addition is performed,
; and the relativistic form is disabled.
; Default: not set (i.e., relativity is applied)
;
; EXAMPLE:
;
; IDL> print, srvadd([0.1d,0,0], [0.5d,0,0])
; 0.56504883 0.0000000 0.0000000
;
; Adds velocities of 0.1 and 0.5 times the speed of light. The
; result is slightly less than the arithmetic sum.
;
;
; IDL> print, srvadd([0.,0.1,0],[0.5d,0,0])
; 0.50000000 0.086602542 0.0000000
;
; Adds velocities in two orthogonal directions. Demonstrates the
; relativistic aberration of velocities (i.e., velocities in the
; perpendicular direction are affected).
;
;
; MODIFICATION HISTORY:
; Written, 28 Jan 2002, CM
; More documentation, 29 Jan 2002, CM
; Add CLASSICAL keyword, 29 Jul 2002, CM
;
; $Id: srvadd.pro,v 1.3 2002/07/29 23:16:47 craigm Exp $
;
;-
; Copyright (C) 2002, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
function srvadd, u, v, classical=classical
nu = n_elements(u)/3
nv = n_elements(v)/3
if nu EQ 0 OR nv EQ 0 then begin
message, 'USAGE: U0 = SRVADD(U1, V)', /info
message, ' U1 = vel. of body in rocket frame; V is velocity of rocket in lab', $
/info
message, ' U0 = vel. of body in lab frame', /info
return, -1d
endif
if nu NE nv AND nu NE 1 AND nv NE 1 then begin
message, 'ERROR: U and V must have the same number of vectors'
endif else if nu EQ 1 then begin
v1 = v
u1 = v*0
u1(0,*) = u(0) & u1(1,*) = u(1) & u1(2,*) = u(2)
endif else if nv EQ 1 then begin
u1 = u
v1 = u*0
v1(0,*) = v(0) & v1(1,*) = v(1) & v1(2,*) = v(2)
endif else begin
u1 = u
v1 = v
endelse
if keyword_set(classical) then begin
return, u1+v1
endif
;; Compute unit vector v, along with 1/gamma and 1 - 1/gamma
vunit = v1
vnorm = total(vunit^2,1) ;; Momentarily = |V|^2
oogamma = sqrt(1 - vnorm) ;; 1/gamma
omoogamma = 1 - oogamma ;; 1 - 1/gamma
vnorm = sqrt(temporary(vnorm)) ;; Now vnorm = |V|
wh = where(vnorm EQ 0, ct) ;; Avoid overflow
if ct GT 0 then vnorm(wh) = 1
;; Normalize the unit vector
if ct GT 0 then $
for i = 0, 2 do vunit(i,*) = vunit(i,*) / vnorm
;; Compute elements of numerator and denominator
udv = total(u1*v1,1) ;; Dot product U1 . V
denom = 1/(1 + udv) ;; Denominator of expression
udvu = temporary(udv)/vnorm ;; Dot product U1 . VUNIT
uu = [1,1,1] ;; Used to expand N-vector to 3xN array
;; U0 = ( (1-1/GAMMA)*(U1 . VUNIT)*VUNIT + U1/GAMMA + V ) / (1 - U1 . V)
return, ((uu#(omoogamma*udvu))*vunit + (uu#oogamma)*u1 + v1)*(uu#denom)
end