We first guess a value for 𝜃, solve for ℎ < (𝑞) [using eq.( )] over a range of wave numbers (we used 𝑞 = {-80000, -79999, … 80000}). We then obtain ℎ(𝑥) from ℎ < (𝑞). With ℎ(𝑥) known, we know tan(𝜃 4O ) and tan(𝜃 4P ) and hence 𝜃 from eq. . We now use this 𝜃 to get a better approximation for ℎ < (𝑞), which will give us better approximation for 𝜃. We see that one gets a stable solution for the shape of our system regardless of initial guess value of 𝜃. The radius corresponding to 𝛬 = 0 is noted as 𝑅 g . For non-zero values of 𝛬, we repeat the above process by varying the radius from 𝑅 g → 𝑅(𝛬) iteratively, till the volume is matched to an accuracy of 10 -6 of relative error.