PRO Dfx_2th, a, esc, dfda, nsc, ma
   
   a1 = [a(0), a(1)]               ;thermal part
   a3 = [a(2), a(3)]               ;thermal part
   f1 = fx_1th(a1, esc)     ;you need the fluxes
   f3 = fx_1th(a3, esc)
   
   dfda = fltarr(nsc, ma)
   
   dfda(*, 0) = f1/a(0)     ;df/de47=f/e47
   dfda(*, 1) = f1*((11.6*esc/a(1)^2)-0.5/a(1))
   ;df/dt6=f*(11.6*e/t6^2-1/(2t6))
   dfda(*, 2) = f3/a(2)     ;df/de47=f/e47
   dfda(*, 3) = f3*((11.6*esc/a(3)^2)-0.5/a(3))
   ;df/dt6=f*(11.6*e/t6^2-1/(2t6))
;also you need dgc/dt6's
   gc1 = 27.83/(a(1)+0.65)^1.33+0.353*a(1)^0.422/esc^0.340
   dg1dt6 = -37.013/(a(1)+0.65)^2.33+0.149/(a(1)^0.578*esc^0.34)
   dfda(*, 1) = dfda(*, 1)+(f1/gc1)*dg1dt6
   
   gc3 = 27.83/(a(3)+0.65)^1.33+0.353*a(3)^0.422/esc^0.340
   dg3dt6 = -37.013/(a(3)+0.65)^2.33+0.149/(a(3)^0.578*esc^0.34)
   dfda(*, 3) = dfda(*, 3)+(f3/gc3)*dg3dt6
   
   RETURN
END