;+
; 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