Specify azimuthal variation

[Zeren Lin]

Hello everyone!

I was wondering if anyone has experience using PyCloudy with (theta, phi) pairs, which means the model is not azimuthal symmetric. I set up a model with theta in the range of (-90, 90) and phi in the range of (0, 360) and run Cloudy for this set pf 1D models. 

I then initialized m3d with no plan_sym and no rotation to test. 

However, the problem was when I tried to plot the emissivities on the 3D cube, the cells in the 3D cube definitely looked like it was cropped: the center was supposed to be (50, 50, 50) but for the z axis this didn't look right.

I looked into the manual and the GitHub code which all said the angle phi is defined as arctan(y/x) + pi but this means it will be slightly different from the (0, 360) convention in spherical coordinates and it has duplicate values. 

Thanks in advance if anyone knows how to solve this problem!

Take care,

Zeren

[Zeren Lin]

Sorry I just noticed that the images were not posted correctly. I attached them here below.

Thanks,

Zeren

[Christophe Morisset]

Hi Zeren,

That's a little bit hard to figure out what is happening without a compete script. Could you generate a simple example and send the python script that leads to the problem you meet?

Thanks a lot,

Christophe 

[Zeren Lin]

Hi Christophe,

Thanks for the swift reply! I will do that and post it later!

Best,

Zeren

[Zeren Lin]

Hi Christophe, 

I've created a Github repository for the project. The .ipynb file is a concise version of the code, the .SED file is the SED profile used, and the .in file is one of the CLOUDY inputs generated by the code. The README has a detailed description of the goal. Basically I was trying to model a galactic disk with the density described by a Sersic-disk profile and the circular velocity set by the halo and disk gravitational potential. I wanted to simulate several line emissions, especially their surface brightness and velocity distribution (and further emission profile). 

I've rerun the code from scratch and noticed that the results are actually symmetric (contrary to what I described in my original post, probably because before I messed up with results from several past runs?). But the problem now is the output emissivities in the 3D cube (as in the final cell) doesn't look like a disk. I was wondering whether this is because the code rescales in the z direction (which is the rotation axis of the disk). 

Please let me know if I've made it clear. Thanks!

Zeren

[Christophe Morisset]

Hi Zeren,

I'm answering in your Github repository in the Issues facility. When the problem is solved, we will post here.

Christophe

[Zeren Lin]

Hi Christophe,

I've made several changes to my code (please take a look at the Simple_rotation_debug4 in Github):

1. The 'negative population and crash' issue was solved by allowing no induced processes in CLOUDY. I'm not totally sure why it is but since I care only about (semi-) forbidden line emissions so I assume it won't make any big differences.

2. I simplified the model a lot realizing I can absorb the axis ratio q later from rotating the grid. Now I only need to specify theta in the range from 0 to 90. I haven't checked whether using the full theta and phi range has the same results but I will if needed later.

3. I generated a plot for the projected velocity and here I have one question: do I understand it correctly that  when defining the velocities in def_profiles_user(), the 3D velocity are specified in the rotated 3D grids instead of the unrotated ones. That is, for my model of a rotating disk viewed from a different angle, I need to rotate my velocities first? (I do have a function for velocity rotation using the same 3D angles in Simple_rotation_debug4).

Thanks!

Zeren

[Christophe Morisset]

Hi Zeren,

I think that this function you wrote to rotate the velocities could be implemented directly in the pyCloudy.c3d.CobCoord object. I'm trying to do it on the 3d_velocities branch of the pyCloudy github repository.

Still in devel, I will advise when completed.

Cheers,

Christophe

[Zeren Lin]

Thanks a lot, Christophe!!

Sorry but I ran into another question while I was trying to interpret the results. Could you kindly help me with that?

In general, I don't quite understand the units of the outputs. For example, when I made a plot as Figure1 attached, I know the units of the XY axes are dimensions of the simulation grids. However, what's the physical units for that? (What should I do to convert it to arcsecs?) Similarly for the Z axis (for the upper plots in Figure1): to generate the plots for the emissivities, I used get_emis function summed along the projection axis (as you did in Example 3), and I was also wondering the units for that? (From the manual, the results from get_emis() is erg/s/cm^3 but it's projected along one axis). If instead I want to get a surface brightness plot to better compare with observations, what should I do?

Thanks,

Zeren

[Christophe Morisset]

Hi Zeren,

In the 3D model, you will have the interpolation of the values as read from the Cloudy outputs. The units are the same. In the case of get_emis, you are right, the unit is erg/s/cm3. If you compute the sum of this cube on a given axis and multiply by the volume of each cell of the 3d cube, you will have erg(s:

m3d.get_emis('O__2_372603A').sum(axis = proj_axis)*m3d.cub_coord.cell_size

This will gives you the intensity emitted by each spaxel of the image. 

If you want to compare with surface brightness, you will need to take into account the distance to the object to transform the geometrical size of the object or the spaxels into angular size.

Christophe

[Christophe Morisset]

Hi Zeren,

I implemented the rotation of the velocity vectors in pyCloudy. I'm not sure it is working well. Could you help in testing that it works well?

You need to update to the branch:

pip install -U git+https://github.com/Morisset/pyCloudy.git@3d_velocities

the set_velocity method should work the same as before. You do not need to rotate the velocity, define it in the "rest frame" without rotation. Then you can set up the angle value and the velocity should turn (I hope...). You can access the (rotated) velocity components with m3d.cub_coord.vel_x (vel_y and vel_z). The norm of the vector is m3d.cub_coord.vel. 

Once it is validated, I will merge this branch and it will be in the master pyCloudy release.

Gimme some feedback,

Cheers,

Christophe

[Zeren Lin]

Hi Christophe,

Sorry for the late reply. I just tested the 3D_velocity branch and it worked as expected!! (Thanks!)

Just to summarize so that others in this group can follow easily:

1. Define the rotation angle

2. Initialize the c3d object with the rotation angle: m3d = pc.C3D(... angles = ...)

3. Define the velocity profile in the rest frame and call it: def_profiles_user(m3d)

This way the velocities are co-rotated with with the 3d cube by the same rotation angle.

============

By the way, as for the units question I had before, do I understand it correctly that since the unit for depths in the Cloudy output files is cm, and I got

m3d.cub_coord.delta_x (same for y and z) = 3.70109031e+21, which is about 1.2 kpc, 

then the units for my plots (even after rotation and projection) are 1.2 kpc for x and y axes?

Thanks a lot,

Zeren

[Christophe Morisset]

Hi Zeren,

Great that it is working! I will merge this branch to the main one.

About the size question, you may check that everything is working as you expect by exploring the Cloudy output files. The #### sequence in the output files is a flag for the first and last zones of the model, where you will find the size of your 1D models. The maximum of the size is what will be used for the 3D mesh. This size divided by the number of cells should give you the size of each cell.

Saludos,

Christophe

[Christophe Morisset]

Hi Zeren,

Well, it seems that the new facility introduced by rotating the velocities does not do the job correctly! The results obtained in the case of a radial expansion with no dependency of the velocity on the angle should not be modified in the new version. This is not the case. Sorry for that, I will have to investigate a little bit more. I reverted the merge and kept the 3d_velocoties branch for testing.

Saludos,

Christophe

[lzc...@gmail.com]

Hi Christophe,

Thanks so much for your help and it has been all really helpful!!

I've managed to make a lot of progress on my simulations for the past months, but I met a few more problems. Recently I've included some other emission lines using emis_tab:

    emis_tab = ['H  1  4861.33A',

            'H  1  6562.81A',

            'Mg 2  2795.53A',

            'Mg 2  2802.71A',

            'O  2  3726.03A',

            'O  2  3728.81A',

            'O  3  5006.84A',

            'Ne 3  3868.76A',

            'Ne 3  3967.47A',

            'Ne 5  3345.99A',

            'Ne 5  3426.03A'

            ]  

Cloudy computed their emissivities successfully but for some reason the 3D model only recognized the [OII] doublet:

m3d.emis_labels

array(['O__2_372603A', 'O__2_372881A'], dtype='<U12')

I believe I'm using the test version you wrote before incorporating velocity rotation, and I was wondering whether you turned off the other elements and only kept [OII] for test purposes. I thought about upgrading to the newest release but I wanted to make sure this is not due to some other reason. I can send you the folder with the simulation results if they are needed. 

PS: in the .in file, there are several 'save last element' commands to save the ionization structure. It seems that they are specified through some config file. Is there any way that I can include Mg there?

save last element hydrogen ".ele_H"

save last element helium ".ele_He"

save last element carbon ".ele_C"

save last element nitrogen ".ele_N"

save last element oxygen ".ele_O"

save last element argon ".ele_Ar"

save last element neon ".ele_Ne"

save last element sulphur ".ele_S"

save last element chlorin ".ele_Cl"

save last element iron ".ele_Fe"

save last element silicon ".ele_Si"

Thanks and take care,

Zeren

[Christophe Morisset]

Hi Zeren,

Well, I'm very sorry, it seems that the full 3D facility development has been moved to old sedimental layers :-( The actual version is not correctly doing the job, as central symmetric models are not behaving like expected. I will try to have a look at the problem.

The way the m3d object works is not the same as the 1D object: the interpolation on the 3D mesh is made only when needed (we do not want to fill the user's memory with a huge amount of data cubes). Then at the beginning there is no emission line cube at all and the emis_label list is empty. On the other side, the emis_labels list of each individual 1D model is filled with all the available labels, you can see it with m3d.m[0].emis_labels which gives you access to the first of your 1D models. 

When you ask for an emission line cube in the 3D model m3d, it is generated from the interpolation and store in the object, for future use: the corresponding line label appears now in the emis_labels list. In your example, you already ask for these 2 OII line cubes, they are stored. The other lines have not been computed, they are not there, but will be computed and stored when asked for.

To add Mg in the element list, change the config variable: pc.config.SAVE_LIST_ELEMS, for example:

pc.config.SAVE_LIST_ELEMS.append(['magnesium', '.ele_Mg'])
print(pc.config.SAVE_LIST_ELEMS)

Tell me if you have other problems,

Cheers,

Christophe

[lzc...@gmail.com]

Hi Christophe,

Thanks for your help! I was able to get emissivities from the other lines successfully!

However, I noticed one other problem (!) when I was writing a proposal and using Cloudy3D to estimate the signal strength... It looks like the outputs such as emissivity, H/e density are properties of the emitting gas clouds. Therefore to calculate properties of the whole nebula like the surface brightness, column densities I need to multiply the local volume filling factor and then sum along the sightline?

Also, you mentioned before the velocity rotation version didn't work for some cases. Could you please explain that a little bit more since it looks like right in my case...

Thanks a lot and stay safe!

Zeren