evaluvating variable inside a function while integrating

edit

To get it more convenient I wrapped numerical_integral for use with piecewise functions:

def numerical_integral_piecewise(f,a,b,**kwds):

    """
    Returns a numerical integral from a piecewise function 

    ``f``  - piecewise function

    ``a``  - lower bound of the interval

    ``b``  - upper bound of the interval

    ``**kwds`` -  options passed to numerical_integral


    EXAMPLES::

        sage: f1(x) = 1
        sage: f2(x) = 2
        sage: f3(x) = 3
        sage: f = piecewise([[(0,4),f1],[(4,6),f2],[(6,10),f3]])
        sage: numerical_integral_piecewise(f,2,9)
        15.0

        if a (or b ) is outside the domain of f, the interval (a,b) is cut 
        to the lower bound (or upper bound) of the domain of f

        sage: numerical_integral_piecewise(f,-2,20)
        20.0

    TODO
        return an error estimate as numerical_integral does

    """    

    D = dict(f.list())
    sum = 0   
    for k in D.keys():
        if k[1] <= a: continue
        if k[0] <= a and b <= k[1]:
            return numerical_integral(D[k],a,b,**kwds)[0]            
        elif k[0] <= a: 
            sum += numerical_integral(D[k],a,k[1],**kwds)[0]
        elif a <= k[0] and k[1] <= b:
            sum += numerical_integral(D[k],k[0],k[1],**kwds)[0]
        elif k[0] <= b <= k[1]:
            sum += numerical_integral(D[k],k[0],b,**kwds)[0]        
    return sum


f1(x) = 1
f2(x) = 2
f3(x) = 3
f = piecewise([[(0,4),f1],[(4,6),f2],[(6,10),f3]])
numerical_integral_piecewise(f,-2,20)

another try:

put the functions in a python dictionary with the respective intervals as keys.

run through the sorted keys and add up the integrals for each fitting interval

from scipy.constants import h, c, k

a=11717

f1(x) = a*(x**-0.53)
f2(x) = a*(x**-0.75)

T = {(1,10000):f1,(10000,20000):f2, (20000,50000):f1}


# Blackbody Planky function
def B(Lambda,Temp):
    return 2*h*c**2/(Lambda**5 *(exp(h*c/(Lambda*k*Temp))-1))

a1 = 300; a2 = 22000  # interval

def flux2(Lambda):
    sum = 0   
    L = copy(T.keys())
    L.sort()
    for k in L:
        if k[1] <= a1: continue
        if k[0] <= a1 and a2 <= k[1]:
            print 'one: key=',k, 'interval from', a1, 'to', a2
            return numerical_integral(2*pi*x*B( Lambda,T[k]),a1,a2)[0]            
        if k[0] <= a1: 
            print 'two: key=',k, 'interval from', a1, 'to', k[1]
            sum += numerical_integral(2*pi*x*B( Lambda,T[k]),a1,k[1])[0]
            continue
        if a1 <= k[0] and k[1] <= a2:
            print 'three: key=',k, 'interval from', k[0], 'to', k[1]
            sum += numerical_integral(2*pi*x*B( Lambda,T[k]),k[0],k[1])[0]
            continue
        if k[0] <= a2 <= k[1]:
            print 'four: key=',k, 'interval from', k[0], 'to', a2
            sum += numerical_integral(2*pi*x*B( Lambda,T[k]),k[0],a2)[0]
        if a1 >= k[1]:
            pass            
    return sum

print 2.5*log(flux2(9000*10**-10))

When calling T2(x) inside the function flux2 x is an element of Symbolic Ring and bool(x < 21500) gives False.

The range 7 to 2150 doesn't matter because T2(x) is evaluated before any substitution of x happens.

You can get the desired behavior with an assumption ( try with and without forget() )

assume(x<=2150)
#print assumptions()
#forget()
from scipy.constants import h, c, k
def T2(x):
    a=11717  
    if x < 21500 : 
        return a*(x**-0.53)
    else :
        print(x.parent())  # Just for debugging        
        return a*(x**-0.75)
# Blackbody Planky function
def B(Lambda,Temp):
    return 2*h*c**2/(Lambda**5 *(exp(h*c/(Lambda*k*Temp))-1))

def flux2(Lambda):
    return numerical_integral(2*pi*x*B(Lambda,T2(x)),7,2150)[0]

print 2.5*log(flux2(9000*10**-10))