Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fitting stream with varying RAR (theta_0, W_0) and disk parameters. #5

Open
martinmestre opened this issue Mar 2, 2022 · 4 comments

Comments

@martinmestre
Copy link
Owner

martinmestre commented Mar 2, 2022

Dear @Charly-Arguelles @danielcarpintero @ankrut :

I quote from the previous issue:

In my opinion, one remaining point would bee to check (at least minimally) what can be gained further regarding the fitting power of the rotation curve, by varying somewhat the baryons (mainly the disk parameters). But of course, it should be re-done the analysis for GD-1, in a kind of iterative fashion (because we are not doing all the fittings in parallel as we have already agreed in the past). Martin, one question from my side (which I forgot): are we using the Miyamoto Nagai disk right?, and, what about the baryonic disks (formula and specific free parameters) used by the other authors we are comparing with in the Vrot curve?. To me, there is certainly some extra freedom about the baryons yet, in the sense that, for example the contribution of the disk could make the RAR DM tail start to decrease further in radius..

I will do some fits varying the barions as follows:
15% for central values of the masses and 10% for characteristic scales, according to the first paragraph of section 3.1.2 of this paper

@martinmestre
Copy link
Owner Author

martinmestre commented Mar 2, 2022

@ankrut @Charly-Arguelles @danielcarpintero , here the results.

The varying parameters in the fit are theta_0, d_theta (=W_0-theta_0) plus the thin and thick disk parameters (Please note that I use the same factors for both disks). The disk parameters are the total mass, the r-scale and the z-scale.

def pot_model(ener_f, theta_0, W_0, beta_0, alpha, scale, scaleb):
    """Potential model including barions and dark matter halo."""
    # Set barionic potentials
    M_gal = 2.32e7  # Msun
    bulge = pot.Plummer(460.0*M_gal, 0.3)
    thin = pot.MiyamotoNagai(1700.0*M_gal*alpha, 5.3*scale, 0.25*scaleb)
    thick = pot.MiyamotoNagai(1700.0*M_gal*alpha, 2.6*scale, 0.8*scaleb)
    # Set halo potential
    param = np.array([ener_f, theta_0, W_0, beta_0])
    halo = pot.RAR(param)
    return bulge, thin, thick, halo

The a priori bounds for those parameters (theta_0, d_theta, alpha, scale, scaleb) are:
# bounds = ((34, 39), (25, 30), (0.85, 1.15), (0.9, 1.1), (0.9, 1.1))
I made two runs with a rather large number of iterations-population equal to 40-40 and 100-100, respectively.
In the first case I obtained:
theta_0, d_theta, alpha, scale, scaleb =
3.642214195068059013e+01,
2.758970214272743249e+01
1.089703158216322576e+00
1.033675113036105708e+00
9.062517097351714401e-01

In the second case I obtained:
3.626724510339099083e+01
2.747475602895485736e+01
1.148997022079204511e+00
1.063196394644405185e+00
9.072897756291734561e-01

Both cases are compatible between them and also in agreement with the case with fixed barions (alpha=1, scale=1, scaleb=1).
The chi^2 reduces a bit to 27 in this two cases with varying disk. In the case of fixed disk, the chi^2 is 29. This slight difference is not noticeable from the stream plots. Please notice that in this case I am using beta_0=1.25e-5 while in the original fits I used 1.2e-5.
Here is a plot of the rotation curve for the second case,
rotation_curve_beta0_1 25e-5
in order to compare with the plots in the previous issue.

@danielcarpintero
Copy link
Collaborator

danielcarpintero commented Mar 2, 2022 via email

@Charly-Arguelles
Copy link
Collaborator

Excellent @martinmestre

@ankrut , @danielcarpintero , as we were discussing today with Martin (from what he figured it out), by defining a NEW likelihood, where we add the \Chi² with the Vrot data (e.g. Sofue+ 2020), to the already used \Chi² of the stream data, it would be very interesting to see the new result of maximization of this new Likelihood, because now we may obtain a considerably better fit to the V_rot , while still providing a good fit to the stream.

As Daniel asks, it would be nice if you can write down how you are planning to define the Lieklihhods (either in the original case, as in the new one)
Carlos

@martinmestre
Copy link
Owner Author

martinmestre commented Mar 2, 2022

@danielcarpintero :

How do you compute chi^2 in these experiments?

Like this:

def chi2(w_0, ener_f, beta_0, ic):
    """Chi^2 function."""
    import wrap
    theta_0 = w_0[0]
    d_theta = w_0[1]
    alpha = w_0[2]
    scale = w_0[3]
    scaleb = w_0[4]
    W_0 = theta_0 + d_theta
    pot_list = pot_model(ener_f, theta_0, W_0, beta_0, alpha, scale, scaleb)
    phi_1, phi_2, d_hel, mu_ra, mu_dec, v_hel, x, y, z, v_circ = orbit_model(
        ic[0], ic[1], ic[2], ic[3], ic[4], ic[5], pot_list)
    cfg.phi_2_spl = interp1d(phi_1, phi_2, kind='cubic')
    cfg.d_hel_spl = interp1d(phi_1, d_hel, kind='cubic')
    cfg.v_hel_spl = interp1d(phi_1, v_hel,  kind='cubic')
    cfg.mu_ra_spl = interp1d(phi_1, mu_ra, kind='cubic')
    cfg.mu_dec_spl = interp1d(phi_1, mu_dec, kind='cubic')
    cfg.phi_1_min = np.amin(phi_1.value)
    cfg.phi_1_max = np.amax(phi_1.value)

    sum = np.zeros(5)

    y_mod = wrap.phi_2_wrap(Iba_sky['phi_1'])
    y_dat = Iba_sky['phi_2']
    sigma2 = 0.5**2
    sum[0] = np.sum((y_dat-y_mod)**2 / sigma2)

    y_mod = wrap.d_hel_wrap(Iba_sky['phi_1'])
    y_dat = Iba_sky['d_hel']
    sigma2 = 1.5**2
    sum[1] = np.sum((y_dat-y_mod)**2 / sigma2)

    y_mod = wrap.mu_ra_wrap(Iba_sky['phi_1'])
    y_dat = Iba_sky['mu_ra']
    sigma2 = 4.0
    sum[2] = np.sum((y_dat-y_mod)**2 / sigma2)

    y_mod = wrap.mu_dec_wrap(Iba_sky['phi_1'])
    y_dat = Iba_sky['mu_dec']
    sigma2 = 4.0
    sum[3] = np.sum((y_dat-y_mod)**2 / sigma2)

    y_mod = wrap.v_hel_wrap(Iba_sky['phi_1'])
    y_dat = Iba_sky['v_hel']
    sigma2 = 100.0
    sum[4] = np.sum((y_dat-y_mod)**2 / sigma2)

    print('chi^2 =', np.sum(sum), ' v_circ = ', v_circ, ' theta_0=', theta_0, ' W_0=', W_0,
          ' alpha=', alpha, ' scale=', scale, 'scaleb=', scaleb)
    return np.sum(sum)

It is a discrete sum of the differences between the Ibata polynomials and our orbits, for the five observables: position in sky (only Phi_2 because Phi_1 is the independent variable here), distance, proper motion and radia velocity, all added together. Likke an L_2 norm between functions.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants