<reponame>bryanmiller/bwm_utils
from __future__ import print_function
import numpy as np

# Python version of the IDL script
# <NAME>

def randomp(pow,n,range_x=[5.,100.]):
# NAME:
#       RANDOMP
# PURPOSE:
#       Generates an array of random numbers distributed as a power law.
# CALLING SEQUENCE:
#       RANDOMP(Pow, N, [ RANGE_X = [low,high]])'
# INPUTS:
#       Pow:  Exponent of power law.
#               The pdf of X is f_X(x) = A*x^pow, low <= x <= high
#               ASTRONOMERS PLEASE NOTE:  
#               pow is little gamma  = big gamma - 1 for stellar IMFs.
#       N:    Number of elements in generated vector.
#
# OPTIONAL INPUT KEYWORD PARAMETER:
#       RANGE_X:  2-element vector [low,high] specifying the range of 
#               output X values# the default is [5, 100].
#
# OUTPUTS:
#       X:    Vector of random numbers, distributed as a power law between
#               specified range
# PROCEDURE:  
#       "Transformation Method" for random variables is described in Bevington 
#       & Robinson, "Data Reduction & Error Analysis for Physical Sciences", 2nd
#       Edition (McGraw-Hill, 1992). p. 83.
#       Output of RANDOMU function is transformed to power-law
#       random variable.
#
# EXAMPLE:
#       Create a stellar initial mass function (IMF) with 10000 stars
#       ranging from 0.5 to 100 solar masses and a Salpeter slope.  Enter:
#
#       MASS = RANDOMP(-2.35,10000,RANGE_X=[0.5,100])
#
# NOTES:
#       Versions 5.1.1 and V5.2 of IDL have a bug in RANDOMU such that the SEED
#       value is initialized to the same value at the start of each session,
#       rather than being initialized by the system clock.    RANDOMP will be
#       affected in a similar manner.
# MODIFICATION HISTORY:
#       Written by <NAME>, Hughes STX, July 13, 1995
#       July 14, 1995   SEED keyword added at Landsman's suggestion.
#                    Documentation converted to standard format.  RSH
#       Converted to IDL V5.0   <NAME>   September 1997
#       Converted to python     <NAME>     October 2014
#-

    if len(range_x) != 2:
        print ('Error - RANGE_X keyword must be a 2 element vector')
        return
        
    pow1 = pow + 1.0
    
    lo = range_x[0]
    hi = range_x[1]
    if lo > hi:
        temp=lo
        lo=hi
        hi=tmp

    r = np.random.random( n )
    if pow != -1.0:
        norm = 1.0/(hi**pow1 - lo**pow1)
        expo = np.log10(r/norm + lo**pow1)/pow1
        x = 10.0**expo
    else:
        norm = 1.0/(np.log(hi) - np.log(lo))
        x = np.exp(r/norm + np.log(lo))

    return x


from __future__ import print_function
import numpy as np

# Python version of the IDL script
# <NAME>

def randomp(pow,n,range_x=[5.,100.]):
# NAME:
#       RANDOMP
# PURPOSE:
#       Generates an array of random numbers distributed as a power law.
# CALLING SEQUENCE:
#       RANDOMP(Pow, N, [ RANGE_X = [low,high]])'
# INPUTS:
#       Pow:  Exponent of power law.
#               The pdf of X is f_X(x) = A*x^pow, low <= x <= high
#               ASTRONOMERS PLEASE NOTE:  
#               pow is little gamma  = big gamma - 1 for stellar IMFs.
#       N:    Number of elements in generated vector.
#
# OPTIONAL INPUT KEYWORD PARAMETER:
#       RANGE_X:  2-element vector [low,high] specifying the range of 
#               output X values# the default is [5, 100].
#
# OUTPUTS:
#       X:    Vector of random numbers, distributed as a power law between
#               specified range
# PROCEDURE:  
#       "Transformation Method" for random variables is described in Bevington 
#       & Robinson, "Data Reduction & Error Analysis for Physical Sciences", 2nd
#       Edition (McGraw-Hill, 1992). p. 83.
#       Output of RANDOMU function is transformed to power-law
#       random variable.
#
# EXAMPLE:
#       Create a stellar initial mass function (IMF) with 10000 stars
#       ranging from 0.5 to 100 solar masses and a Salpeter slope.  Enter:
#
#       MASS = RANDOMP(-2.35,10000,RANGE_X=[0.5,100])
#
# NOTES:
#       Versions 5.1.1 and V5.2 of IDL have a bug in RANDOMU such that the SEED
#       value is initialized to the same value at the start of each session,
#       rather than being initialized by the system clock.    RANDOMP will be
#       affected in a similar manner.
# MODIFICATION HISTORY:
#       Written by <NAME>, Hughes STX, July 13, 1995
#       July 14, 1995   SEED keyword added at Landsman's suggestion.
#                    Documentation converted to standard format.  RSH
#       Converted to IDL V5.0   <NAME>   September 1997
#       Converted to python     <NAME>     October 2014
#-

    if len(range_x) != 2:
        print ('Error - RANGE_X keyword must be a 2 element vector')
        return
        
    pow1 = pow + 1.0
    
    lo = range_x[0]
    hi = range_x[1]
    if lo > hi:
        temp=lo
        lo=hi
        hi=tmp

    r = np.random.random( n )
    if pow != -1.0:
        norm = 1.0/(hi**pow1 - lo**pow1)
        expo = np.log10(r/norm + lo**pow1)/pow1
        x = 10.0**expo
    else:
        norm = 1.0/(np.log(hi) - np.log(lo))
        x = np.exp(r/norm + np.log(lo))

    return x
# Python version of the IDL script # <NAME> # NAME: #       RANDOMP # PURPOSE: #       Generates an array of random numbers distributed as a power law. # CALLING SEQUENCE: #       RANDOMP(Pow, N, [ RANGE_X = [low,high]])' # INPUTS: #       Pow:  Exponent of power law. #               The pdf of X is f_X(x) = A*x^pow, low <= x <= high #               ASTRONOMERS PLEASE NOTE: #               pow is little gamma  = big gamma - 1 for stellar IMFs. #       N:    Number of elements in generated vector. # # OPTIONAL INPUT KEYWORD PARAMETER: #       RANGE_X:  2-element vector [low,high] specifying the range of #               output X values# the default is [5, 100]. # # OUTPUTS: #       X:    Vector of random numbers, distributed as a power law between #               specified range # PROCEDURE: #       "Transformation Method" for random variables is described in Bevington #       & Robinson, "Data Reduction & Error Analysis for Physical Sciences", 2nd #       Edition (McGraw-Hill, 1992). p. 83. #       Output of RANDOMU function is transformed to power-law #       random variable. # # EXAMPLE: #       Create a stellar initial mass function (IMF) with 10000 stars #       ranging from 0.5 to 100 solar masses and a Salpeter slope.  Enter: # #       MASS = RANDOMP(-2.35,10000,RANGE_X=[0.5,100]) # # NOTES: #       Versions 5.1.1 and V5.2 of IDL have a bug in RANDOMU such that the SEED #       value is initialized to the same value at the start of each session, #       rather than being initialized by the system clock.    RANDOMP will be #       affected in a similar manner. # MODIFICATION HISTORY: #       Written by <NAME>, Hughes STX, July 13, 1995 #       July 14, 1995   SEED keyword added at Landsman's suggestion. #                    Documentation converted to standard format.  RSH #       Converted to IDL V5.0   <NAME>   September 1997 #       Converted to python     <NAME>     October 2014 #-