Cook and Felderman analysis to perform data reduction on heated surfaces This is a "pain", but I wonder if I can get any real speedup with the follow Python script - it is taking well over 10 minutes to calculate the way it is written now. (And that's only when feeding around 1-2 promille of the data).

I have an with shape (roughly): [260, 300, 300] 

Now I am doing a cougar=feldermann analysis on this:

import numpy as np
import math
#COOK AND FELDERMAN
def CandF(dT, time, rho, c, k):
    sqrt = np.sqrt
    sumdata=np.zeros(len(time))
    for n in range(len(time)):    
        for data in range(1,n+1):
            sumdata[n] += dT[data] / (sqrt(time[n]-time[data])+sqrt(time[n]-time[data-1]))

    qs = (2*sqrt(rho*c*k))/(sqrt(math.pi)) * sumdata
    return qs

def data_reduction_constant_time(dT, timestep, rho, cp, k):
    #print (len(dT), timestep)
    time = np.linspace(0, len(dT)*timestep, int(len(dT)))
    return data_reduction(dT, time, rho, cp, k )

def data_reduction(dT, time, rho, cp, k):
    return CandF(dT, time, rho, cp, k)

class data:
    def __init__(self):
        #normally lots of file reading

        self.integration = 400
        self.cols = 320
        self.rows = 260
        self.depth = 300
        self.readData()

    def readData(self):
        self.ml_delta_temp = 3 * np.random.random_sample((self.cols, self.rows, self.depth))

    def calcQML(self):
        print("--- Calculating q ---")
        print (self.cols, self.rows, self.depth)
        v = np.zeros([self.cols,self.rows,self.depth])
        T = 590
        P = 28 * 100000
        R = 287
        c = 1005
        k = 1.4
        q = 20
        rho = P/(R*T)  
        print(v)

        print("-------------")
        for x in range(self.cols ):
            print ("calculating column: ", x)
            for y in range(self.rows ):
                v_tmp = data_reduction_constant_time(self.ml_delta_temp[x,y,:], 
                                                                self.integration,
                                                                rho, c, k)
                v[x,y,:] = v_tmp
        print("--- q done ---")
        print(v)

if (__name__ == "__main__"):
    dat = data()
    dat.calcQML()


self.ml_delta_temp is a (260, 300, 300) shaped numpy array. cols, and rows are the size of the array, and integration is the step size (how much time is between each 3rd axis slice).

Now the problem is: on my computer it takes almost a minute to calculate a single row.

How can I speed it up, or as the same question: in what part is this function losing the most time?

I can show the mathematical formula for the Cook and Felderman equation:

$$qCF = \frac{2\sqrt{\rho ck}}{\sqrt{\pi}}\sum_{i=1}^n \frac{T_i-T_{i-1}}{\sqrt{t_n-t_i}+\sqrt{t_n-t_{i+1}}}$$

Berry, S.A. and Horvath, T.J. Discrete
Roughness Transition for Hypersonic Flight Vehicles. Journal of Spacecraft and Rockets, 2:216{
227, March-April 2008. DOI 10.2514/1.30970.

Now this is to get the total heat at time "t_n", I am however interested in the heat flux profile, how it changes over time. So I need to calculate this formula for each n. Also I do know not T, but only delta-T (but as one can see that is only beneficial).
I suggest the following replacement for function CandF():

def CandF3(dT, time, fac):
sqrt = np.sqrt
sumdata = np.zeros(len(time))
for n in range(len(time)):
    tn = time[n]
    # precompute sqrt(tn-time[0])
    last_sq = sqrt(tn - time[0])
    sum_n = 0
    for j in range(1, n + 1):
        sq_j = sqrt(tn - time[j])
        sum_n += dT[j] / (sq_j + last_sq)
        last_sq = sq_j
    sumdata[n] = sum_n
return sumdata * fac


The most rewarding idea here was that from the 2 square roots per innermost loop one is the second is computed repeatedly and can be saved for the next loop (here: last_sq). Halving the number of np.sqrt() calls nearly halves the runtime for large numbers of time steps (len(time)).
More optimizations are the use of a local variable for a loop-invariant (tn), use of a scalar when indexed variable is not indexed incrementally (sum_n), and precomputing the final constant sqrt(rho*c*k*2.0) in the caller. This seems natural as the contributing constants are all defined there.