2017-07-19 20 views
1

私はいくつかのデータを持っており、私は波長(青の点)の大きさをプロットしました。モデルの星雲をファイルから読み込み、同じグラフ(ピンクの線)にプロットするコードを作成します。このコードでは、この線をグラフ上で上下に動かすように調整できるスケールがあります。これまで私は、私のポイントに目を向けるほど線が近くなるように尺度を変更していましたが、尺度の値を計算するコードを書きたいと思います。ラインは最小です。このように見えるグラフを与える点と曲線の間の最小距離を見つけるためのPythonコード

#Import modules 

from math import * 
import numpy as np 
import matplotlib.pyplot as plt 

# Specify data 

wavelength = 
np.array([357.389,445.832,472.355,547.783,620.246,752.243,891.252,2164.089]) 
magnitude = 
np.array([24.0394,23.1925,23.1642,22.4794,21.7496,20.9047,20.4671,19.427]) 

# Create Graph 

#plt.scatter(wavelength, magnitude) 
#plt.ylim([25,18]) 
#plt.xlim([300,2200]) 
#plt.xlabel('wavelength (nm)') 
#plt.ylabel('magnitude') 
#plt.title('object 1') 
#plt.show() 
#plt.close() 

#now - here is some code that reads a model stellar population model from a 
file 

lines = open('fig7b.dat').readlines() 

wavelengths, luminosities = [],[] 

for l in lines: 
    s = l.split() 
    wl = s[0] 
    old = s[-1] 
    if '#' not in wl: 
     wavelengths.append(float(wl)) #wavelength in angstroms 
     luminosities.append(float(old)) #luminosities are in log units! 


scale = 3.5 
c=3.e8 
wavelengths = np.array(wavelengths) 
nus = c/(wavelengths*1.e-10) 
luminosities = np.array(luminosities) + scale 

luminosity_density = np.log10(((10**luminosities)*wavelengths)/nus) 

#plt.plot(wavelengths,luminosity_density) 
#z = 1.0 
#plt.plot(wavelengths*(1+z),luminosity_density,color='r') 

#plt.axis([900, 10000, 25,31]) 
#plt.savefig('sed.png') 
#plt.show() 
#plt.close() 

Mpc_to_cm = 3.086e24 #convert Mpc to cm 
z = 0.3448 #our chosen redshift 
D_L = 1841.7 * Mpc_to_cm 

#remember luminosity_density is logged at the moment 
flux_density = (10**luminosity_density) * (1+z)/(4*pi*D_L**2) #units will 
be erg/s/cm^2/Hz 

#now turn that into an AB magnitude - goes back to log 
AB_mag = -2.5*np.log10(flux_density) - 48.6 

#try plotting your photometry on here and play with z and D_L 
plt.plot(wavelengths*(1+z),AB_mag,color='pink') 
plt.scatter(wavelength*10., magnitude,color='cornflowerblue') 
plt.axis([900, 25000, 30,18]) 
plt.xlabel('wavelength') 
plt.ylabel('magnitude') 
plt.title('object 1') 
plt.savefig('sed_ab.png') 
plt.show() 

enter image description here

はまた、最高のスケール値を印刷すると便利だろうこれは、これまでの私のコードです。 私は非常にPythonとプログラミング一般的に新しく、ピンクの線は単純な方程式ではありません(私が与えられたファイルには多くのデータ点が含まれています)ので、ちょっと立ち往生しています。私が問題を記述するのに正しい言語を使用していない場合や長いコードについては、謝罪しています。多くのコメントは、私が別のプロットをしていたときに私の上司が以前から守っていたものです。 (私のpython 2.7を使用しています)

リンクfig7b.datする:https://drive.google.com/open?id=0B_tOncLLEAYsbG8wcHJMYVowOXc

+0

モデル曲線のデータ点のRMSを計算し、['scipy.optimize.minimize'](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize)で最小化することができます.minimize.html)。また、[フィット感](https://en.wikipedia.org/wiki/Goodness_of_fit)もご覧ください。 –

+0

'fig7b.dat'のコピーを受け取る可能性はありますか? –

+0

@HughBothwell私は一番下のリンクをアップロードしました! –

答えて

1

まず、各点は点の最初のリストに対応するように(それぞれの対応する対を曲線データから点のリストを作成します点は同じX座標、すなわち同じ波長を持つ)。

これらの2つのポイントの最小距離は、単に(sum(points2)-sum(points1))/len(points1)です。このウィルプリント値2.25を実行

次の例を見て

points1 = [1.1, 1.4, 1.8, 1.9, 2.3, 1.7, 1.9, 2.7] 
points2 = [8.4, 3.5, 2.9, 7.6, 0.1, 2.2, 3.3, 4.8] 

def min_distance(first,second): 
    assert len(first) == len(second) # must have same size 
    result = (sum(second) - sum(first))/len(first) 
    return result 

print("Adding this value to the first series of points") 
print("will provice minimum distance between curves") 
print(min_distance(points1,points2)) 

。 をpoints1のすべての値に追加すると、2つのポイントのセット間の最小距離(この場合は62.36)が得られます。

問題では、points1magnitudeの配列になります。 points2は、波長に対応するfig7b.datからのポイントになります。

これは、点と曲線の間の平方和を最小にしたいと仮定しています。また、距離が垂直に測定されていると仮定しています(そのため、対応する波長の点を抽出する必要があります)。あなたがspicy.optimizeを使用せず、独自の小さなコードを書きたい場合は

1

私が推薦する :

は、あなたの観測波長のそれぞれの理論値を評価するために、あなたの理論スペクトルの補間を使用しますhttps://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html

例えば:あなたは\カイを計算することができるよりも

from scipy.interpolate import interp1d  
f2 = interp1d(wavelengths, luminosities, kind='cubic') 

^あなたが試したいと思っているすべてのスケール値に対して{2}を計算し、その後に最小値を見つけます。

関連する問題