Small Breakthrough in the Flux code

After all the bureaucracy I had to solve after coming back to Brasil, and a small vacation of one week in the middle of nowhere, last week the work restarted for good...

I had a small problem trying to recompile my codes here in the Astrophysics Lab, but it is not worth mentioning here. The Fortran compiler here was a different one from the compiler I was using in Ann Arbor. After some small changes in the makefile, all was good again.

I am still following the guidelines I described in the last post here. So, after compiling the Flux codes again, I decided to add the wind optical depth component to the stellar continuum. After the calculation was done, I had a unpleasant surprise. All the absorption due to the disk wind was red shifted instead of being blue shifted. That was totally wrong and unexpected. Actually, the amount of absorption due to the wind in the red side of the spectrum should be almost non-existent. After that I was sure there was a problem in the velocity field projected onto the line of sight.

The routine that calculates the velocity projections is called 'chi.f', and is also the routine that calculates the opacity and emissivity of each point in the system. Here are the lines originally used by this routine to calculate the velocity projections:

...
vzd = vz
cosphi = xd/sqrt(xd**2+yd**2)
sinphi = yd/sqrt(xd**2+yd**2)
vyd = vr*sinphi+vphi*cosphi

if (vyd.eq.0) then
arg = teta
else
arg = pi/2 - teta - arctan(vzd/vyd)
endif

vel = sqrt(vyd**2+vzd**2)

if (yd.lt.0) vel = -1 * vel

vlos = vel*cos(arg)
...

The first thing I noticed about these lines was a small inconsistency that does not change actually the results. It is the fact that these lines first project the velocities in the y-axis, instead of the x-axis. The light rays are always parallel to the x-axis, so the velocities should also be projected in the x-axis. But we also have azimuthal symmetry, so it really does not matter in which direction in xy-plane the velocity is projected. I corrected it anyway, just for consistency. So this line

vyd = vr*sinphi + vphi*cosphi,

was modified to

vxd = vr*cosphi - vphi*sinphi.

I started to play with these projections a little bit, and the plot the distribution of points with vlos greater than zero, or lesser than zero, in a 3D-graph, and then a projection of these points in xy-plane. I noticed that some changes in the signs in the projection equations, gave me the blue absorption component that I was looking for, but in all the cases the point distributions where not as I would expect. After making some play, and some tests, I decided to follow each of those lines in the code, starting with vzd, plotting the 3D distribution of points with vzd over zero, and bellow zero. The same for vr, then vphi, the vxd, and so on...

All those distributions were correct, but vlos was wrong. Then I saw the line that was wreaking all the havok...

if (yd.lt.0) vel = -1 * vel

The code was using the modulus of vel, that is always positive. But the projected velocity is not always positive, so there should be a change of signs in the velocity somewhere. First, I have changed the projection from the y-axis, to the x-axis, so instead of yd lesser than zero, the change of signs should happen for xd lesser than zero. But that was not the only problem... We have to consider rotation in the system, so the axis where there will be the change of velocities is also rotated. In the case of no rotation, simply changing yd for xd in the above mentioned line would do the trick, but in the case with rotation that won't work...

To calculate the velocity onto the line of sight, the code calculates the poloidal velocity first (vel), then it changes the sign of vel, and then it projects onto the line of sight. To solve the sign problem, I bypassed the poloidal velocity calculation as it is really not needed anywhere else in the code, and then projected vzd and vxd directly onto the line of sight. I just saved a small amount of time bypassing some unneeded lines of code, simplifying the code, and making the projection directly the problem with the velocity sign disappears. The modified code is now here:

...
vzd = vz
cosphi = xd/sqrt(xd**2+yd**2)
sinphi = yd/sqrt(xd**2+yd**2)
vxd = vr*cosphi - vphi*sinphi

vlos = vzd*cos(teta)+vxd*sin(teta)
...

All the other lines were commented, and bypassed. After these modifications the distribution of points with vlos greater than zero, or lesser than zero, seems to be correct, and now the wind absorption feature is in the blue side of the profile, as expected.