function randomp,n,indx=indx,range=range,seed=seed

; ----------------------------------------------------------
;+
; NAME:
;       RANDOMP
;
; PURPOSE:
;       Generate array of random deviates with power law distribution
;
; AUTHOR:
;       Simon Vaughan (U.Leicester)
;
; CALLING SEQUENCE:
;       y = randomp(100,indx=-1.5,range=[0.1,1.0],seed=seed)
;
; INPUTS:
;       n     - (integer) length of output array
;
; OPTIONAL INPUTS:
;       indx  - (float) power law index
;       range - (vector) [min,max] values for output deviates
;       seed  - (integer) seed for random number generator
;               identical results on different runs.
; OUTPUTS:
;       y     - (vector) array of n random numbers
;
; DETAILS:  
;       Generates an array of N random deviates that follow a
;       power law probability density function (PDF) with an
;       index INDX between values RANGE[0] and RANGE[1].
;       If no values are given for N, INDX or RANGE then default
;       ones are used.
;
;       The function uses the "transformation method" - also known
;       as "inverse transform sampling" - to convert random deviates
;       distributed uniformly on [0,1] to a power law distribution.
;       See "Numerical Recipes" by W. H. Press et al. section 7.2.
;
; SOLUTION:
;       We aim to transform from uniform x to power law y.
;
;         p(x) = 1           0 < x < 1
;         p(y) = N*y^b      y0 < y < y1
;
;       From calculus the transformation from x to y is: 
;
;         |p(y)dy| = |p(x)dx|
;         p(y) = |dx/dy| p(x)
;
;       We must therefore solve dx/dy = f(y). The solution is x = F(y)
;       and its inverse gives y = F^-1(x) which is the transformation
;       we seek. From the above we know that:
;
;         p(y)dy = |dx/dy| dy (because p(x)=1 over 0 < x < 1) 
;
;       and for a power law:
;
;         N*y^b = dx/dy = f(x)
;
;       If we integrate this (w.r.t. y) we get x = f(y)dy = F(y)
;
;         x = F(y) = N/(b+1)*y^(b+1) + C          for b /= -1
;         x = F(y) = N*ln(y) + C                  for b =  -1
;
;       where C is the constant of integration.
;       Re-arrange this to get y = F^-1(x):
;
;         y = F^-1(x) = [(x-C)*(b+1)/N]^[1/(b+1)] for b /= -1
;                     = exp[(x-C)/N]              for b =  -1
;
;       We must find N and C from the boundary conditions.
;       Firstly, we can find N because the integral of p(x)dx
;       from x=0 to x=1 must be 1. This gives 
;         
;         N = (b+1)/[y_max^(b+1) - y_min^(b+1)]   for b /= -1
;         N = 1/[ln(y_max) - ln(y_min)]           for b =  -1
;
;       And we also know that y_min = F^-1(x=0) so therefore:
;
;         C = N*y_min^(b+1)                       for b /= -1
;         C = N*ln[y_min]                         for b =  -1
; 
;       And finally, substituting for N and C we get:
;
;         y = [x*(b+1)/N + y_min^(b+1)]^[1/(b+1)] for b /= -1
;         y = exp[x/N + ln y_min]                 for b =  -1
;
;       These are the two equations used for the transformation
;       from uniform to power law distributed random deviates
;       (depending on whether the power law index is -1 or not).
;
; HISTORY:
;       23/07/2007  -- v1.0 -- first working version
;-
; ----------------------------------------------------------

; watch out for errors
  on_error,2

; ----------------------------------------------------------
; Check the arguments and input

; if N is not set use default
  if (n_elements(n) eq 0) then n=1

; if INDX not set use default
  if (n_elements(indx) eq 0) then indx=0

; if RANGE not set use default
  if (n_elements(range) eq 0) then range=[0.01,1.0]

; check RANGE makes sense
  if (range[0] eq range[1]) then begin
      print,'** RANDOMP: MIN/MAX values are the same in RANGE keyword'
      return,0
  endif

; ----------------------------------------------------------
; Main section

; set lower/upper limits on Y
  if (range[0] lt range[1]) then begin
      y0 = range[0]
      y1 = range[1]
  endif else begin
      y0 = range[1]
      y1 = range[0]
  endelse     

; drawn N uniform random deviates
  x = randomu(seed,n)

; perform transformation from X -> Y where X is
; uniformly distributed as Y is power law distributed.
; INDX=-1 is a special case that is handled differently.
; First calculate in log-space, then convert to linear.

  if (indx ne -1.0) then begin
      n0 = 1.0/(y1^(indx+1) - y0^(indx+1))
      logy = 1.0/(indx+1.0) * alog10(x/n0 + y0^(indx+1))
      y = 10.0^logy
  endif else begin
      n0 = 1.0/(alog(y1) - alog(y0))
      y = exp(x/n0 + alog(y0))
  endelse

; ----------------------------------------------------------
; Return the result to the user
  return, y

END