import numpy as np import scipy as sp import matplotlib.pyplot as plt plt.interactive(True) ''' function returns one of 4 additive contributions to v2, square of the rotation curve of a galaxy. Explicit arguments are the radius (array), scale lengths of bulge, disk and halo. Mass of a component implicitly equal to 1e6 Msun, via appropriate scaling of gravitational constant G. ''' def v2(n,r,rb,rd,rh): if(n==0): # central superm. black hole M=1e6 Msun v2 = G/r return(v2) if(n==1): # Plummer sphere bulge, M_b=1e6 Msun, core radius rb v2 = G*r*r/(r*r+rb*rb)**1.5 return(v2) if(n==2): # disk of mass=Msun & expon. scale rd, Md=1e6 Msun v2 = G/r*(1-(1+r/rd)*np.exp(-r/rd)) return(v2) if(n==3): # dark halo, Mc = 4 pi r_h^3 rho_h = 1e6 Msun v2 = G/rh - (G/r)*np.arctan(r/rh) return(v2) else: print('***err') return(0*r-1) ''' total v2(r) array assembled according to mass array m[0:4] ''' def v2_tot(r,rb,rd,rh): sum = 0*r for p in range(N): sum += m[p]*v2(p,r,rb,rd,rh) return(sum) ''' Given a set of scale lengths r_bulge, r_disk and r_halo, this function solves for 4 paramaters (masses in units of 1e6M_sun), which multiply the appropriate unit v2(r) contributions. The function also returns the best fit v_fit array and value E of the best chi-square, a measure of quadratic misfit. ''' def find_m (r_b,r_d,r_h): # fill A and b for n in range(N): v2n = v2(n,r,r_b,r_d,r_h) b[n] = np.sum(v2n*YY) # create b_n for l in range(n+1): # array A[l,n] dotprod = np.sum(v2(l,r,r_b,r_d,r_h)*v2n) A[l,n] = dotprod A[n,l] = dotprod # solve NxN matrix equation A a = b m = np.linalg.solve(A,b) # best fit v_fit = v2_tot(r,r_b,r_d,r_h)**0.5 # misfit (chi-square, really) E = np.mean(((v_fit-Y)/sig)**2) return(m,E) ''' Pyplot displays color image of arr with title and labels ''' def show_arr(arr, title, xlab, ylab, clean=0): if(clean): plt.cla() x = np.linspace(0.9*rh,2.1*rh,np.size(arr,axis=1)) y = np.linspace(0.9*rd,1.6*rd,np.size(arr,axis=0)) plt.figure(figsize=(10.5,8.5)) # or: dpi=140 # rescale array c = plt.pcolor(x,y,arr) b = plt.colorbar(c,orientation='vertical') # q = plt.winter() # a_max = np.maximum(arr) e = plt.contour(x,y,arr,np.linspace(1.,1.30,30)**4,colors='k') #plt.axis(xlabel=xlab,ylabel=ylab) plt.xlabel(xlab, fontsize=18) plt.ylabel(ylab, fontsize=18) plt.title(title, fontsize=18) plt.text(rh_best,rd_best,'+', fontsize=18) plt.text(8.4,4.72,'best fit r$_d$='+str(np.around(rd_best,3)),fontsize=14) plt.text(15,2.75,'best fit r$_h$='+str(np.around(rh_best,3)),fontsize=14) plt.show() input(' continue? ') ''' Main program fit-glx-4+2.py Fit multi-parameter galaxy rotation curve model (4+2 parameters) to synthetic observations of rotatian curve of a disk galaxy, consisting of 4 components: (i) supermassive black hole (BH) of mass m[0] in the center, (ii) central bulge of unknown core radius r_b= and core mass m[1] (modeled as Plummer sphere), (iii) exponential-surface density disk of unknown scale length r_d and total mass m[2], and (iv) dark matter halo modeled as pseudo-isothermal sphere of core radius r_h and core mass m[3]. M = 80 = number of synthetic observational data points Unit of mass is 1e6 M_sun Unit of length is kpc Units of rotational velocity are km/s Known: r_b = 1 kpc Unknowns: r_d, r_h, and 4 masses M_BH, M_bulge, M_disk, M_halo (core masses) placed in array m[]. Method: The only embedded parameters, which cannot be found directly by Gauss's Least Squares fit because they are not front coefficients of additive functions, are r_d and r_h. We have to search for them by some other algorithm, whereas the N masses (for a fixed set of r_b,r_d,r_h) are quickly found by an NxN system of linear algebraic equations. This hybrid method avoids an expensive 6-parameter non-linear search, splitting the search into a 4-parameter calculation by linear algebra, and a slower 2-parameter search on top of 4-parameter calculation. ''' plt.figure(dpi=140) N = 4 # number of unknown masses M = 80 # number of observations rb, rd, rh = 1, 3, 9 # length scales of components of synthetic model G = 4.302 # kpc/(10^6 Msun)*(km/s)^2, value sig = 3. # km/s, obs. uncert. A = np.zeros((N,N),dtype=float) m = np.zeros(N) b = np.zeros(N) r = np.linspace(0.1,30,M) # hidden truth about scale lengths rr = [rb,rd,rh] # hidden truth about masses M_hidden = np.array([5e2,1e4,1e5,1.5e5]) # M:BH,b,d,h [10^6 Msun] m = M_hidden v2_waves = 60**2 *np.cos(r*(45-r)/45)*(r/15)*np.exp(-r/15) # some waves v_hidden = (v2_tot(r,rb,rd,rh)+v2_waves)**0.5 # synthetic model # square of observational data Y = v_hidden +sig*np.random.randn(M) # hidden model + noise YY = Y*Y plt.scatter(r,Y,s=8,color=(0,0,1),linewidth=1) plt.grid() plt.xlabel('r [kpc]') plt.ylabel('V [km/s]') plt.title('4+2 param. model of galaxy as SMBH+bulge+disk+DM_halo') plt.show() input(' ok?') plt.plot(r,(v_hidden**2-v2_waves)**0.5,color=(0.5,0.4,0.8),linewidth=1) plt.plot(r,v_hidden,color=(1,0,1),linewidth=0.9) plt.show() input(' ok?') err_min = 1e9 #print('hidd. rh %8.4f m %7.1f %7.0f %7.0f %7.0f E %8.5f' # %(rh, m[0],m[1],m[2],m[3], E_find)) K = 100 chi2arr = np.zeros((K,K),dtype=float) # double loop fills KxK array of chi2's and finds minimum chi2 for i in range(K): r_d = rd*0.9 + i*rd*0.7/K for j in range(K): r_h = rh*0.9+ j*rh*1.2/K # for fixed disk and halo radii r_d & r_h, find 4 masses m[] and chi2 m, chi2 = find_m (rb,r_d,r_h) chi2arr[i,j] = chi2 # store all chi2 for plotting if (err_min > chi2): rd_best = r_d; i_best = i rh_best = r_h; j_best = j err_min = chi2; m_best = m # what's wrong with he first column j=0? print('chi2[0:2,0:K/4]\n',chi2arr[0:2,0:K//4],i_best,j_best) # print found parameters print('best: %3d %3d rd%7.4f rh %7.4f m:%7.1f %7.0f %7.0f %7.0f chi2%8.5f' %(i_best,j_best,rd_best,rh_best,m[0],m[1],m[2],m[3], err_min)) # build rotation curve for the best fit m = m_best v_best = v2_tot(r,rb,rd_best,rh_best)**0.5 # plot total & separate contributions to the best rotation curve v(r) for i in range(N): plt.plot(r,(m[i]*v2(i,r,rb,rd_best,rh_best))**0.5, color=[i/3,.6,1-i/3],linewidth=1.5) plt.plot(r,v_best,color=(1,0,0),alpha=0.6,linewidth=1.5) plt.grid() fs = 10 plt.text(13.7,220,'rotation curve',fontsize=fs) plt.text(14,125,'halo',fontsize=fs) plt.text(14,111,str(np.around(m[3]/1000,1))+'$\cdot 10^9 M_\odot$',fontsize=fs) plt.text(7,190,'disk',fontsize=fs) plt.text(7,176,str(np.around(m[2]/1000,1))+'$\cdot 10^9 M_\odot$',fontsize=fs) plt.text(23,66,'bulge',fontsize=fs) plt.text(23,52,str(np.around(m[1]/1000,2))+'$\cdot 10^9 M_\odot$',fontsize=fs) plt.text(15,15,'black hole',fontsize=fs) plt.text(21,15,str(np.around(m[0]/1000,3))+'$\cdot 10^9 M_\odot$',fontsize=fs) plt.plot([0,30],[0, 0],color=(.6,.6,.6),alpha=0.6,linewidth=0.9) plt.plot([0,0],[0,260],color=(.6,.6,.6),alpha=0.6,linewidth=0.9) plt.show() #input(' overplot smooth assumed model components?') #plt.plot(r,v_hidden,color=(0,.75,1),alpha=0.5,linewidth=1) #for i in range(N): # plt.plot(r,(M_hidden[i]*v2(i,r,rb,rd,rh))**0.5, # color=[i/6,.2,1/2-i/6],linewidth=0.7) #plt.grid() #plt.show() input(' done? chi2: ') # plot error field print(chi2arr[20:40,20:40]) chi2arr = chi2arr/np.min(chi2arr) chi2arr = np.minimum(chi2arr,2.0) txt = '$\chi^2$ map (min = '+str(np.around(err_min,3))+')' show_arr (chi2arr,txt,'r$_h$','r$_d$') # compute velocity dispersion in dark halo # sigma^2 = G M_h / r_h by virial theorem sigma_h = np.around((G *m[3]/r_h)**0.5, 3) print(' velocity dispersion of dark halo core is',sigma_h,' km/s') # #sp.optimize.minimize(fun, x0, ..)