I am trying to have contourplots of deviation dc for variables a and b.

```
In(27):= (Alpha)1(l_, b_, a_, (Theta)_, M_,
r_) := -(a^2*(1 + l)*(2*M + 2*r + b) +
r*(2*r^2 + 3*b*r + b^2 - 2*M*(3*r + b)))/((a*
Sqrt(1 + l)*(2*M - 2*r - b))*Sin((Theta)))
```

```
In(28):= (Beta)1(l_, b_, a_, (Theta)_, M_, r_) :=
Sqrt((r^2*(8*a^2*(1 + l)*
M*(2*r + b) - (2*r^2 + 3*b*r + b^2 - 2*M*(3*r + b))^2)/((1 +
l)*a^2*(2*M - 2*r - b)^2)) +
a^2*Cos((Theta))^2*(1 + l) - (Alpha)1(l, b, a, (Theta), M, r)^2*
Cos((Theta))^2)
```

```
In(29):= re1(l_, b_, a_, (Theta)_, M_) :=
RankedMax(
r /. NSolve({(Beta)1(l, b, a, (Theta), M, r) == 0, r > 0}, r,
Reals), 1)
```

```
In(30):= pr1(l_, b_, a_, (Theta)_, M_) :=
RankedMax(
r /. NSolve({(Beta)1(l, b, a, (Theta), M, r) == 0, r > 0}, r,
Reals), 2)
```

```
In(31):= xc1(l_, b_, a_, (Theta)_, M_) :=
NIntegrate((Alpha)1(l, b, a, (Theta), M, r)*(Beta)1(l, b,
a, (Theta), M, r)*D((Alpha)1(l, b, a, (Theta), M, r), r), {r,
pr1(l, b, a, (Theta), M), re1(l, b, a, (Theta), M)})/
NIntegrate((Beta)1(l, b, a, (Theta), M, r)*
D((Alpha)1(l, b, a, (Theta), M, r), r), {r,
pr1(l, b, a, (Theta), M), re1(l, b, a, (Theta), M)})
```

```
In(32):= R(l_, b_, a_, (Theta)_, M_,
r_) := ((Alpha)1(l, b, a, (Theta), M, r) -
xc1(l, b, a, (Theta), M))*
D((Beta)1(l, b, a, (Theta), M, r), r) - (Beta)1(l, b,
a, (Theta), M, r)*D((Alpha)1(l, b, a, (Theta), M, r), r)
```

```
In(33):= raverage(l_, b_, a_, (Theta)_, M_) :=
Sqrt((1/(Pi))*
NIntegrate(
R(l, b, a, (Theta), M, r), {r, re1(l, b, a, (Theta), M),
pr1(l, b, a, (Theta), M)}))
```

```
In(34):= R1(l_, b_, a_, (Theta)_, M_,
r_) := ((Sqrt(((Alpha)1(l, b, a, (Theta), M, r) -
xc1(l, b, a, (Theta), M))^2 + (Beta)1(l, b, a, (Theta),
M, r)^2) -
raverage(l, b, a, (Theta),
M))^2)*((((Alpha)1(l, b, a, (Theta), M, r) -
xc1(l, b, a, (Theta), M))*
D((Beta)1(l, b, a, (Theta), M, r), r) - (Beta)1(l, b,
a, (Theta), M, r)*
D((Alpha)1(l, b, a, (Theta), M, r),
r))/(((Alpha)1(l, b, a, (Theta), M, r) -
xc1(l, b, a, (Theta), M))^2 + (Beta)1(l, b, a, (Theta), M,
r)^2))
```

```
In(41):= dc1(l_, b_, a_, (Theta)_,
M_) := (1/raverage(l, b, a, (Theta), M))*
Sqrt((1/(Pi))*
NIntegrate(
R1(l, b, a, (Theta), M, r), {r, re1(l, b, a, (Theta), M),
pr1(l, b, a, (Theta), M)}));
```

```
In(43):= ContourPlot(
ConditionalExpression(dc1(0, b, a, 17*(Pi)/180, 1),
Im(dc1(0, b, a, 17*(Pi)/180, 1)) == 0), {b, 0, .42^2}, {a, .001, 1})
```

```
Out(43)= $Aborted
```

The code calculates the average radius of a parametric curve and then find the deviation of the curve from circularity in terms of rms value. the last line is not giving any output neither giving any error.