It is been a year since I last wrote here. Since two weeks ago, the flux code is finally working, and flawlessly. I spent a lot of time trying to make things to work as I would like them to do, and it took me more time than I wanted and expected, but now it is working the way I expect.
The thing that took most of my time solving was the multi-level atom. The method I was using to calculate the 3-level hydrogen atom was flawed, and I was using it just to make the code run. Now I am using the output source function value of CV to calculate a 3-level hydrogen atom. CV uses a 2-level atom with continuum approach to calculate the line source function, without needing the number densities of the required levels, but using the line source function is possible to acquire the number densities. Since the medium is considered to be optically thin for Lyman lines, the population of the fundamental level can be calculated using thermodynamical equilibrium. Knowing the first level number density, and the ratio of number densities between the line transition levels, we can solve a 3-level atom. It seems simple, but this took months of my time. I tried many different ways to calculating these number densities, before I saw the way I am using now.
After that was done, I had to solve some problems with the code that were causing problems with the results. There were problem when calculating the grid used by Flux. Problems with the size of the grid boundaries. A problem after the grid re-sampling, when points that should be out of the region of diskwind and magnetospheric-accretion have values of densities, velocities, etc, different than zero. And the last problem solved, was that one of the codes was replacing the phi-component of the velocity to zero everywhere.
Now the code is running, the grid points used are correct, and also is the velocity and density fields. No changes were made in the temperature law.
Tuesday, May 5, 2009
Tuesday, April 8, 2008
Stellar continuum plus diskwind flux
After solving the problem with the velocity projections onto the line of sight, the next step was to add the full diskwind contribution to the stellar flux. I had to solve some small problems in the code, but nothing worth of note here. After some debugging and coding, I got my first line profile calculated with the 'Flux' code that contains the stellar continuum and the diskwind contribution.
Comparing this plot with the last post's plot, we can see that the wind contribution is far superior than the stellar continuum, but still we have a blue-shifted absorption feature, that is expected.
Comparing this plot with the last post's plot, we can see that the wind contribution is far superior than the stellar continuum, but still we have a blue-shifted absorption feature, that is expected.
Thursday, April 3, 2008
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.
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.
...
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.
Thursday, March 6, 2008
Stellar Continuum with Flux
I will try another approach to solve the flux calculation problem I am having in the codes. I will calculate each component of the profile separately, doing one step at a time. First only the continuum coming from the star, then see what happens with the continuum under the optical depths of the magnetosphere, and the wind, but without the flux contribution from these sources. In this way I can see if the absorptions on the continuum are in the right place, and the velocities are right. Then I will calculate only the magnetospheric contribution to the profile, and I will do the same to the diskwind component. Then I can see where the broad wings are coming from. After figuring that out, I will finally and both component to the stellar continuum, and I hope to have a good profile at last.
During my first continuum calculation I noticed a problem that was happening for velocities just higher than zero velocity. I forgot to bypass the optical depth test, and with that the code bypassed all the continuum calculations from those velocities. After ignoring the optical depth test in 'intens.f', I have the right continuum again. I guess this problem might be linked to the always present deep red shifted absorption (more than expected, and every time) found in all the profiles calculated so far. I need to take a deeper look at that.
The continuum calculated is in the file named 'prof.60.ha.mmax_7700K.dmax_7700K.continuum'. That is for a 60 degree inclination, the h-alpha line, maximum temperature in the magnetospheric funnel of 7700K, and maximum temperature in the diskwind region of 7700K. The stellar parameters are for a typical Classica T Tauri Star:
Obs: Just for you to remember. The intens.f have been changed to calculate only the continuum. I guess you will not forget that, but anyway...
During my first continuum calculation I noticed a problem that was happening for velocities just higher than zero velocity. I forgot to bypass the optical depth test, and with that the code bypassed all the continuum calculations from those velocities. After ignoring the optical depth test in 'intens.f', I have the right continuum again. I guess this problem might be linked to the always present deep red shifted absorption (more than expected, and every time) found in all the profiles calculated so far. I need to take a deeper look at that.
The continuum calculated is in the file named 'prof.60.ha.mmax_7700K.dmax_7700K.continuum'. That is for a 60 degree inclination, the h-alpha line, maximum temperature in the magnetospheric funnel of 7700K, and maximum temperature in the diskwind region of 7700K. The stellar parameters are for a typical Classica T Tauri Star:
- Mass: 0.5 Msun
- Radius: 2.0 Rsun
- Photospheric Temperature: 4000K
- Accretion Ring Temperatura: 8000K
- Mass Accretion Rate: 1e-8 Msun/Year
- Mass Loss Rate: 1e-9 Msun/Year
- Magnetospheric Inner and Outer Radius: 3.0 Rstar and 3.3 Rstar
- Diskwind Inner and Outer Radius: 3.31 Rstar and 20.0 Rstar
Obs: Just for you to remember. The intens.f have been changed to calculate only the continuum. I guess you will not forget that, but anyway...
Wednesday, February 20, 2008
Fortran and Fluxes Problems
In the last few weeks I have been working in this "Flux" code non-stop. I have been solving problems trying to make that code work with the output of my modified CV code for a little more than two months, and I still was not able to make it work properly.
Last week I found a problem, that I first thought was an interpolation problem, but it turned out to be something else. Actually, the problem was something that I have done wrong. When I tested the points through each of the rays to find if they were inside the wind region, I made a mistake. The wind has some geometric parameters, but to use those parameters I first have to made a change of variables, and I was not doing that, so my wind region was being calculated with wrong boundaries. And these wrong boundaries were the reason I had problems with the interpolation.
I recalculated the correct boundaries, the interpolation problems were solved, but still the code did not work properly.
After a while I noticed other thing wrong. My grid that should be a circular grid, was instead an elliptical grid. I solved that today. The grid is circular, but the results are far from the expected yet.
I am running out of ideas to solve the problems with this code... I have the feeling that my stellar flux is far too low compared with the magnetosphere and diskwind fluxes, but I don't know why, and that does not make much sense to me. The optical depths in the magnetospheric funnel seems to be right this time. But there is a weird red-shifted absorption feature in my profiles that does not match what I was expecting for the inclination I am using.
May the Radiative Tranfer and Fortran 77 Gods have mercy on my mind, and give me some inspiration to solve this... :P
Last week I found a problem, that I first thought was an interpolation problem, but it turned out to be something else. Actually, the problem was something that I have done wrong. When I tested the points through each of the rays to find if they were inside the wind region, I made a mistake. The wind has some geometric parameters, but to use those parameters I first have to made a change of variables, and I was not doing that, so my wind region was being calculated with wrong boundaries. And these wrong boundaries were the reason I had problems with the interpolation.
I recalculated the correct boundaries, the interpolation problems were solved, but still the code did not work properly.
After a while I noticed other thing wrong. My grid that should be a circular grid, was instead an elliptical grid. I solved that today. The grid is circular, but the results are far from the expected yet.
I am running out of ideas to solve the problems with this code... I have the feeling that my stellar flux is far too low compared with the magnetosphere and diskwind fluxes, but I don't know why, and that does not make much sense to me. The optical depths in the magnetospheric funnel seems to be right this time. But there is a weird red-shifted absorption feature in my profiles that does not match what I was expecting for the inclination I am using.
May the Radiative Tranfer and Fortran 77 Gods have mercy on my mind, and give me some inspiration to solve this... :P
Subscribe to:
Posts (Atom)