title: Plot 3D isosurface in python
body: <p>I am trying to plot the following function in python using plotly or matplotlib for a given value of omega:</p>
  
  <pre><code>omega = (1/(12*np.pi*r**3)*((3*np.cos(THETA)**2-1)+1.5*np.sin(THETA)**2*np.cos(2*PHI)))
  </code></pre>
  
  <p>To do this I specify the value of <code>omega</code>, calculate <code>r</code> and then convert from polar to cartesian coordinates.</p>
  
  <pre><code>import numpy as np
  import plotly.graph_objects as go
  
  omega = 10
  theta, phi = np.linspace(0,2*np.pi, 400), np.linspace(0,np.pi, 400)
  THETA,PHI = np.meshgrid(theta,phi)
  
  #Calculate R for a given value of omega
  R1 = (1/(12*np.pi*omega)*((3*np.cos(THETA)**2-1)+1.5**np.sin(THETA)**2*np.cos(2*PHI)))
  R = np.sign(R1) * (np.abs(R1))**(1/3)
  
  #Convert to cartesians
  X = R * np.sin(PHI) * np.cos(THETA)
  Y = R * np.sin(PHI) * np.sin(THETA)
  Z = R * np.cos(PHI)
  
  #Plot isosurface plot
  fig = go.Figure(data=[go.Surface(z=Z, x=X, y=Y)])
  fig.show()
  </code></pre>
  
  <p>However, by doing this I lose the information about the positive and negative lobes, i.e. here both are plotted.</p>
  
  <p><a href="https://i.stack.imgur.com/A8Xiu.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/A8Xiu.png" alt="enter image description here"></a></p>
  
  <p>I can get round this by setting negative values of omega to <code>NaN</code>. This switches off the negative lobes but results in rendering artefacts for the graphs.</p>
  
  <pre><code>import numpy as np
  import plotly.graph_objects as go
  
  omega = 10
  theta, phi = np.linspace(0,2*np.pi, 400), np.linspace(0,np.pi, 400)
  THETA,PHI = np.meshgrid(theta,phi)
  
  #Calculate R for a given value of omega
  R1 = (1/(12*np.pi*omega)*((3*np.cos(THETA)**2-1)+1.5**np.sin(THETA)**2*np.cos(2*PHI)))
  R = np.sign(R1) * (np.abs(R1))**(1/3)
  
  #Remove negative lobes
  R[R1 &lt; 0.] = np.NaN
  
  #Convert to cartesians
  X = R * np.sin(PHI) * np.cos(THETA)
  Y = R * np.sin(PHI) * np.sin(THETA)
  Z = R * np.cos(PHI)
  
  #Plot isosurface plot
  fig = go.Figure(data=[go.Surface(z=Z, x=X, y=Y)])
  fig.show()
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/1Y3rZ.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/1Y3rZ.png" alt="enter image description here"></a></p>
  
  <p>I'm not sure how to overcome this issue - if I increase the number of points for <code>THETA</code> and <code>PHI</code>, the graph renders very slowly, and it's still not possible to increase the number of points sufficiently to remove the artefacts completely. Ideally I would pass <code>r</code>, <code>theta</code> and <code>phi</code> values to the plotting function and plot an isosurface at a given value of <code>omega</code>, but this is possible only in cartesians. Converting the whole function to cartesians would lead to <code>f(x,y,z,r</code>) which I also can't plot.</p>