mirror_int_fac = exp(vpa^2 * (mudB/dz)/(mudB/dz + Zedpihnc/dz)) is the integrating factor needed to turn the dg/dvpa part of the GKE advance into an advection equation a, b and c contain the sub-, main- and super-diagonal terms, respectively if running in full-flux-surface mode, solve mirror advance in y-space rather than ky-space due to alpha-dependence of coefficients corresponds to sign of mirror term positive on RHS of equation must treat boundary carefully assumes fully upwinded at outgoing boundary corresponds to sign of mirror term negative on RHS of equation must treat boundary carefully assumes fully upwinded at outgoing boundary time_upwind = 0.0 corresponds to centered in time time_upwind = 1.0 corresponds to fully implicit (upwinded) account for fact that we have expanded d(gnorm)/dvpa, where gnorm = g/exp(-v^s); this gives rise to d(gnormexp(-vpa^2))/dvpa + 2vpagnormexp(-vpa^2) term we solve for gnorm*exp(-vpa^2) and later multiply by exp(vpa^2) to get gnorm multiply by mirror coefficient