[Theta] = .24;
    sourcex = 5.57;
    energysource = 1.35;
    energysoft = .48;
    energyhard = 2.8;
    iterations = 12;
    tg = Table[
       Arg[Conjugate[
           FSg[j, [Theta], sourcex, energysoft, energysource,
            iterations]]*
          FHg[j, [Theta], sourcex, energyhard, energysource,
           iterations]]/wRs[[j]], {j, 1, 12}];
    freq = {.000075, .00015, .00025, .0003, .0007, .001, .0015, .0022,
    .0032, .0045, .007, .01};
    comp = Thread[{freq, tg}];
    source = {{.000075, 65}, {.00015, 165}, {.00025, 125}, {.0003,
        35}, {.0007, -10}, {.001, -20}, {.0015, -30}, {.0022, -15},
    {.0032, -20}, {.0045, -5}, {.007, -3}, {.01, -8}};
    ListLogLinearPlot[{comp, source}, PlotRange -> {{10^(-5), .02}, All},
     PlotStyle -> {Red, Black}, InterpolationOrder -> 1, Frame -> True,
     PlotLegend -> {"theory", "observations"},
     LegendPosition -> {-.7, -.3}, LegendSize -> {.6, .2},
     PlotMarkers -> Automatic, LegendShadow -> None,
     FrameLabel -> {"Temporal Frequency (Hz)", "Lag(s)"},
     ]]]

For[j = 1, j < 12, j++,
 FSg[j_, theta_, r0_, energys_, energy0_, nmax_] :=
  Module[{stuffs, zetamins, zetamaxs, finalvals, zetas, zeta0, sums,
    terms, [Mu], [Kappa]},
   zetas = energys/511.25/theta; zeta0 = energy0/511.25/theta;
   zetamins = Min[zetas, zeta0];
   zetamaxs = Max[zetas, zeta0]; [Kappa] =
    1/2*((3*khat*vhat)/(2*theta) + 4);
   stuffs = (-3*N0*khat*Exp[(I*wtilda[[j]]*y0)]*
       Exp[((zeta0 - zetas)/2)]*zetas^([Kappa] - 4))/(2*[Pi]*
       theta^4*(8.19*10^(-7))^3*(r0)^(1/2)*(Rs)^3*zeta0^([Kappa]));
   sums = 0;
   Do[
    [Mu] =
     1/2*((3 - (3*khat*vhat)/(2*theta))^2 +
         4*(3*khat*vhat)/(2*theta)*(2*sdal[[j, n]])/vhat)^.5;
    terms =
     Gamma[[Mu] - [Kappa] + 1/2]/(IMH[[j, n]]*Gamma[1 + 2*[Mu]])*
      listhOUT[j, n, r0]*listhOUT[j, n, NRs]*
      WhittakerM[[Kappa], [Mu], zetamins]*
      WhittakerW[[Kappa], [Mu], zetamaxs];
    sums = sums + terms,
    {n, 1, nmax}];
   finalvals = sums*stuffs; finalvals
   ]]

For[j = 1, j < 13, j++,
 FHg[j_, theta_, r0_, energyh_, energy0_, nmax_] :=
  Module[{stuffh, zetaminh, zetamaxh, finalvalh, zetah, zeta0, sumh,
    termh, [Mu], [Kappa]},
   zetah = energyh/511.25/theta; zeta0 = energy0/511.25/theta;
   zetaminh = Min[zetah, zeta0];
   zetamaxh = Max[zetah, zeta0]; [Kappa] =
    1/2*((3*khat*vhat)/(2*theta) + 4);
   stuffh = (-3*N0*khat*Exp[(I*wtilda[[j]]*y0)]*
       Exp[((zeta0 - zetah)/2)]*zetah^([Kappa] - 4))/(2*[Pi]*
       theta^4*(8.19*10^(-7))^3*(r0)^(1/2)*(Rs)^3*zeta0^([Kappa]));
   sumh = 0;
   Do[
    [Mu] =
     1/2*((3 - (3*khat*vhat)/(2*theta))^2 +
         4*(3*khat*vhat)/(2*theta)*(2*sdal[[j, n]])/vhat)^.5;
    termh =
     Gamma[[Mu] - [Kappa] + 1/2]/(IMH[[j, n]]*Gamma[1 + 2*[Mu]])*
      listhOUT[j, n, r0]*listhOUT[j, n, NRs]*
      WhittakerM[[Kappa], [Mu], zetaminh]*
      WhittakerW[[Kappa], [Mu], zetamaxh];
    sumh = sumh + termh,
    {n, 1, nmax}];
   finalvalh = sumh*stuffh; finalvalh
   ]]