PRO Dfx_thppl, a, esc, dfda, nsc, ma
   
   a3 = [a(0), a(1)]               ;thermal part
   a1 = [a(2), a(3)]               ;non-thermal part
   f3 = fx_1th(a3, esc)     ;you need the fluxes
   f1 = fx_1pl(a1, esc)
   
   dfda = fltarr(nsc, ma)
   
   dfda(*, 0) = f3/a(0)     ;df/e47=f/e47
   dfda(*, 1) = f3*((11.6*esc/a(1)^2)-0.5/a(1))
;df/dt6=f*(11.6*e/t6^2-1/(2t6))
   dfda(*, 2) = f1/a(2)     ;df/dk=f/k
   dfda(*, 3) = -f1*alog(esc)   ;df/dgamma=-ln(e)*f
   
;also you need dgc/dt6
   gc = 27.83/(a(1)+0.65)^1.33+0.353*a(1)^0.422/esc^0.340
   dgdt6 = -37.013/(a(1)+0.65)^2.33+0.149/(a(1)^0.578*esc^0.34)
   dfda(*, 1) = dfda(*, 1)+(f3/gc)*dgdt6
   
   RETURN
END