function ISCO_eq_plus(x::Float64,f::Field, g::Field, h::Field, W::Field,rh::Float64,static=false)
    #ISCO EQUATION TO OBTAIN LOCAITON OF CO-ROTATING ORBIT
    y=pi/2
    if static
        return (2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gphiphi(x,y,f,g,h,W,rh) - 2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh) - gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh) + gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))/(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gphiphi(x,y,f,g,h,W,rh) - gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh))
    else
        return (-2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^4*gphiphi(x,y,f,g,h,W,rh)^2*gtphi(x,y,f,g,h,W,rh) - 8*gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^5*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)^2 + 8*gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^4*gtphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)*(2*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gphiphi(x,y,f,g,h,W,rh) + gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh)) + gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + 2*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*(-(gtphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)) + gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh))*sqrt((gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2 - gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true))*(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh) - gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh))^2*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))^2) + 2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*(-(gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh)) + gtphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtt(x,y,f,g,h,W,rh))*sqrt((gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2 - gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true))*(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh) - gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh))^2*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))^2) + gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtphi(x,y,f,g,h,W,rh)*(-2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh)^2 + (gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh) - 4*gtphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh) + gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtt(x,y,f,g,h,W,rh))*(-gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))) + 2*gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^3*gphiphi(x,y,f,g,h,W,rh)^2*gtphi(x,y,f,g,h,W,rh) - 3*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)^2 + 2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtphi(x,y,f,g,h,W,rh)*(2*gtphi(x,y,f,g,h,W,rh)^2 - 5*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh)*(gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh) + gtphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtt(x,y,f,g,h,W,rh))*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))*(gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtphi(x,y,f,g,h,W,rh) + gtphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh) + 2*gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))) + gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^3*(gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtphi(x,y,f,g,h,W,rh)*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*(-8*gtphi(x,y,f,g,h,W,rh)^3 + 4*gphiphi(x,y,f,g,h,W,rh)*gtphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))) - 2*gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^3*(2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh)*(4*gtphi(x,y,f,g,h,W,rh)^2 - 3*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gphiphi(x,y,f,g,h,W,rh)*(4*gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + gtt(x,y,f,g,h,W,rh)*(gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + gtt(x,y,f,g,h,W,rh)*(gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh) - gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)))) + gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*(-2*gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtt(x,y,f,g,h,W,rh)*sqrt((gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2 - gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true))*(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh) - gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh))^2*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))^2) + 2*gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*sqrt((gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2 - gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true))*(-(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh)) + gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh))^2*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))^2) + gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh)^2*(-gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + 2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^3*gphiphi(x,y,f,g,h,W,rh)*(4*gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^2*(-4*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh)*(-4*gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + (-gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))*(2*gtphi(x,y,f,g,h,W,rh)*(gtphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh) + gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh)) + gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))) + gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*(2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^3*gtt(x,y,f,g,h,W,rh)^3 - gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*(-gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))*(-2*gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh)^2 + gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh) - 6*gtphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh) + gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtt(x,y,f,g,h,W,rh)^2)))) / ((-(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh)) + gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh))*(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gphiphi(x,y,f,g,h,W,rh)^2 + gtt(x,y,f,g,h,W,rh)*(4*gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gphiphi(x,y,f,g,h,W,rh) - 4*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh) + gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh)) - 2*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*(2*gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtphi(x,y,f,g,h,W,rh) - 2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))))
    end
end

function ISCO_eq_minus(x::Float64,f::Field, g::Field, h::Field, W::Field,rh::Float64,static=false)
    #ISCO EQUATION TO OBTAIN LOCAITON OF COUNTER-ROTATING ORBIT
    y=pi/2
    if static
        return (2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gphiphi(x,y,f,g,h,W,rh) - 2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh) - gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh) + gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))/(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gphiphi(x,y,f,g,h,W,rh) - gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh))
    else
        return (-2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^4*gphiphi(x,y,f,g,h,W,rh)^2*gtphi(x,y,f,g,h,W,rh) - 8*gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^5*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)^2 + 8*gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^4*gtphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)*(2*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gphiphi(x,y,f,g,h,W,rh) + gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh)) + 2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*(gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh) - gtphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtt(x,y,f,g,h,W,rh))*sqrt((gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2 - gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true))*(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh) - gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh))^2*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))^2) - gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*(2*(-(gtphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)) + gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh))*sqrt((gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2 - gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true))*(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh) - gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh))^2*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))^2) + gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)*(-gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))) + gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtphi(x,y,f,g,h,W,rh)*(-2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh)^2 + (gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh) - 4*gtphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh) + gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtt(x,y,f,g,h,W,rh))*(-gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))) + 2*gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^3*gphiphi(x,y,f,g,h,W,rh)^2*gtphi(x,y,f,g,h,W,rh) - 3*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)^2 + 2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtphi(x,y,f,g,h,W,rh)*(2*gtphi(x,y,f,g,h,W,rh)^2 - 5*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh)*(gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh) + gtphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtt(x,y,f,g,h,W,rh))*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))*(gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtphi(x,y,f,g,h,W,rh) + gtphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh) + 2*gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))) + gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^3*(gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtphi(x,y,f,g,h,W,rh)*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*(-8*gtphi(x,y,f,g,h,W,rh)^3 + 4*gphiphi(x,y,f,g,h,W,rh)*gtphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))) - 2*gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^3*(2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh)*(4*gtphi(x,y,f,g,h,W,rh)^2 - 3*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gphiphi(x,y,f,g,h,W,rh)*(4*gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + gtt(x,y,f,g,h,W,rh)*(gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + gtt(x,y,f,g,h,W,rh)*(gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh) - gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)))) + gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*(2*gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtt(x,y,f,g,h,W,rh)*sqrt((gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2 - gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true))*(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh) - gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh))^2*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))^2) - 2*gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*sqrt((gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2 - gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true))*(-(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh)) + gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh))^2*(gtphi(x,y,f,g,h,W,rh)^2 - gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))^2) + gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh)^2*(-gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + 2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^3*gphiphi(x,y,f,g,h,W,rh)*(4*gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^2*(-4*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh)*(-4*gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh)) + (-gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))*(2*gtphi(x,y,f,g,h,W,rh)*(gtphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh) + gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh)) + gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))) + gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*(2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^3*gtt(x,y,f,g,h,W,rh)^3 - gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*(-gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))*(-2*gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh)^2 + gtt(x,y,f,g,h,W,rh,dx=2,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh) - 6*gtphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh) + gphiphi(x,y,f,g,h,W,rh,dx=2,dr=true)*gtt(x,y,f,g,h,W,rh)^2))))/((-(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh)) + gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtt(x,y,f,g,h,W,rh))*(gtt(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gphiphi(x,y,f,g,h,W,rh)^2 + gtt(x,y,f,g,h,W,rh)*(4*gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gphiphi(x,y,f,g,h,W,rh) - 4*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh) + gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)^2*gtt(x,y,f,g,h,W,rh)) - 2*gtt(x,y,f,g,h,W,rh,dx=1,dr=true)*(2*gtphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtphi(x,y,f,g,h,W,rh) - 2*gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gtphi(x,y,f,g,h,W,rh)^2 + gphiphi(x,y,f,g,h,W,rh,dx=1,dr=true)*gphiphi(x,y,f,g,h,W,rh)*gtt(x,y,f,g,h,W,rh))))
    end
end

#function ISCO(f::Field, g::Field, h::Field, W::Field,rh::Float64,q::Float64=0.0; guess_minus = nothing, guess_plus = nothing)
function ISCO(f::Field, g::Field, h::Field, W::Field,rh::Float64,q::Float64=0.0; guess = nothing)
    y=pi/2
    M=GetMass(f,g,h,W,rh)
    j=GetJ(f,g,h,W,rh)
    chi = j/M^2
    rhK=M/2*sqrt(1-chi^2-q^2)

    static=false
    if chi < 1e-10
        static=true
    end

    #SEE https://duetosymmetry.com/tool/kerr-isco-calculator/
    z1=1+(1-chi^2)^(1/3) * ( (1+chi)^(1/3) + (1-chi)^(1/3) )
    z2=sqrt(3*chi^2+z1^2)
    rbl_minus = (3+z2+sqrt((3-z1)*(3+z1+2*z2)))
    rbl_plus= (3+z2-sqrt((3-z1)*(3+z1+2*z2)))

    rbl_minus=M*find_zero(r->4*(q^2-r)*(sqrt(r-q^2)+chi)^2+r*(q^2+(r-2)*r+chi^2),rbl_minus)
    rbl_plus=M*find_zero(r->4*(q^2-r)*(sqrt(r-q^2)-chi)^2+r*(q^2+(r-2)*r+chi^2),rbl_plus)

    xkerr_plus=rBLKerr2x(rbl_plus,rhK,chi,q)
    xkerr_minus=rBLKerr2x(rbl_minus,rhK,chi,q)

    #=
    if isnothing(guess_minus)
        guess_minus = xkerr_minus
        #guess_minus = rBLKerr2x(rbl_minus/2.0,rhK,chi,q)
    end
    if isnothing(guess_plus)
        guess_plus = xkerr_plus
        #guess_plus = rBLKerr2x(rbl_plus*8.0,rhK,chi,q)
    end
    =#
    if isnothing(guess)
        guess_minus = xkerr_minus
        guess_plus = xkerr_plus
    else
        guess_plus = guess[1]
        guess_minus = guess[2]
    end

    sol_plus = fzero(x->ISCO_eq_plus(x,f,g,h,W,rh,static),guess_plus)
    sol_minus = fzero(x->ISCO_eq_minus(x,f,g,h,W,rh,static),guess_minus)

    #FOR KERR WE HAVE M*w+- = +- 1/((rBL+-/M)^(3/2) +- χ). SCHWARZSCHILD = 1/(6*√6)
    w_plus = (-gtphi(sol_plus,y,f,g,h,W,rh,dx=1,dr=true)+sqrt(gtphi(sol_plus,y,f,g,h,W,rh,dx=1,dr=true)^2-gtt(sol_plus,y,f,g,h,W,rh,dx=1,dr=true)*gphiphi(sol_plus,y,f,g,h,W,rh,dx=1,dr=true)))/gphiphi(sol_plus,y,f,g,h,W,rh,dx=1,dr=true)
    w_minus = (-gtphi(sol_minus,y,f,g,h,W,rh,dx=1,dr=true)-sqrt(gtphi(sol_minus,y,f,g,h,W,rh,dx=1,dr=true)^2-gtt(sol_minus,y,f,g,h,W,rh,dx=1,dr=true)*gphiphi(sol_minus,y,f,g,h,W,rh,dx=1,dr=true)))/gphiphi(sol_minus,y,f,g,h,W,rh,dx=1,dr=true)

    #wKerr_plus = 1/((rbl_plus/M)^(3/2)+chi)/M
    #wKerr_minus = -1/((rbl_minus/M)^(3/2)-chi)/M

    wKerr_plus= (rbl_plus^2*sqrt(M*(rbl_plus-M*q^2)) + M^3*q^2*chi - M^2*rbl_plus*chi) / (rbl_plus^4+M^3*(M*q^2-rbl_plus)*chi^2)
    wKerr_minus= (-rbl_minus^2*sqrt(M*(rbl_minus-M*q^2)) + M^3*q^2*chi - M^2*rbl_minus*chi) / (rbl_minus^4+M^3*(M*q^2-rbl_minus)*chi^2)

    RKerr_plus = 2*rhK*sqrt(gKerrN(xkerr_plus,y,rhK,0.0,chi,q)/((-1+xkerr_plus)^2*fKerrN(xkerr_plus,y,rhK,0.0,chi,q)))
    RKerr_minus = 2*rhK*sqrt(gKerrN(xkerr_minus,y,rhK,0.0,chi,q)/((-1+xkerr_minus)^2*fKerrN(xkerr_minus,y,rhK,0.0,chi,q)))

    R_plus = CircumferencialRadius(sol_plus,f,g,h,W,rh)
    R_minus = CircumferencialRadius(sol_minus,f,g,h,W,rh)

    println()
    println("ISCO DETAILS")
    println("Prograde (co-rotating) Circular Orbit of Massive Particles Located at x = ", sol_plus, ". r/rh = ", 2/(1-sol_plus), ". Circumferencial Radius/M = ", R_plus/M, ". Coordinate shift w.r.t. KerrN (x-xKerrN) = ", sol_plus-guess_plus, ". Difference to Comparable KerrN (R/RKerrN-1) = ", R_plus/RKerr_plus-1,  ". ω*M = ", w_plus*M, ". ω/ωKerrN - 1 = ", w_plus/wKerr_plus-1)
    println("Retrogade (counter-rotating) Circular Orbit of Massive Particles Located at x = ", sol_minus, ". r/rh = ", 2/(1-sol_minus), ". Circumferencial Radius/M = ", R_minus/M, ". Coordinate shift w.r.t. KerrN (x-xKerrN) = ", sol_minus-guess_minus, ". Difference to Comparable KerrN (R/RKerrN-1) = ", R_minus/RKerr_minus-1, ". ω*M = ", w_minus*M, ". ω/ωKerrN - 1 = ", w_minus/wKerr_minus-1)
    println()

    return [R_plus/RKerr_plus-1, w_plus/wKerr_plus-1, R_minus/RKerr_minus-1, w_minus/wKerr_minus-1, sol_plus, sol_minus]
end