Sharley法によるバックグラウンドのプログラム

前提・実現したいこと

Shirley(シャ―リ法）で複数ピークを持つような波形（Si2pなど）にバックグラウンドを引きたい。

発生している問題・エラーメッセージ

計算した値をlist内に入れていきたいのですが、listの中に値が入らないので、そもそも自作したプログラム自体が間違っているのではないかと思っています。

該当のソースコード

Python
1
2import numpy as np
3import  matplotlib.pyplot as plt
4
5x,y = np.loadtxt ("Si2p.txt",skiprows=1,unpack = True)
6
7#shirley法によるバックグラウンドのプログラム
8Blist = []
9y1 = np.sum(y) - y[-1]
10h = 0.05
11
12
13k = 1
14Flist =[]
15while True:
16    for i in range(203):
17        R1 = h * (y1 - 0.5*(y[1] - y[-1]) + (y[k] - y[-1]))
18        Q1 = h * (y1 - 0.5*(y[k] - y[-1] + y[203] - y[-1]))
19        B1 = (y[1] - y[203])*(Q1/R1) + y[-1]
20        F = y[k] - B1
21        k += 1
22        Blist.append(B1)
23        Flist.append(F)
24    for i in range(100):
25        t = 1
26        z = Flist[t]
27        Q = h * (y1 - 0.5*(z - y[-1] + y[203] - y[-1]))
28        R = h * (y1 - 0.5*(y[1] - y[-1]) + (z - y[-1]))
29        B = (y[1] - y[203])*(Q/R) + y[-1]
30        Blist.append(B)
31        Flist.append(F)
32        t += 1                
33             
34print(Blist)
35
36plt.scatter(x,y,marker="o" ,color="black",s=10)
37ax = plt.gca()        
38y_min, y_max = ax.get_ylim()
39ax.set_ylim(0,6000)
40x_min, x_max = ax.get_xlim()
41#ax.set_xlim(97,102)
42plt.title("Shirley",fontname="MS Gothic")
43plt.plot(x,y,color="black")
44plt.xlabel("BandEnergy[eV]",fontname="MS Gothic")
45plt.ylabel("Intensity[a.u.]",fontname="MS Gothic")
46ax.invert_xaxis()
47
48plt.show()
49

試したこと

プログラムを複数個自分の力で一から書いて試してみたが、どれも同じようにlistに値を入れて更新するいう作業がうまくいかない。Pandasで読み込んでいたのをnumpyのみにかえるなど様々な事を試してみたがうまくいかないので、プログラム自体が間違っている可能性が高い。

補足情報（FW/ツールのバージョンなど）

以下のリンクにあるシャーリー法の計算式をプログラムに起こそうと奮闘しています。一か月ぐらい試してやっているのですが、時間だけ溶けてしまうので、もしプログラムに詳しい方がいればアドバイスお願いします。
初めての投稿でおかしな部分があるとは思いますがご容赦ください。
https://www.jstage.jst.go.jp/article/jvsj2/56/6/56_13-LC-019/_pdf
[Si2p data]
x y
106.95 100
106.9 200
106.85 329
106.8 315
106.75 323
106.7 309
106.65 333
106.6 307
106.55 315
106.5 311
106.45 290
106.4 317
106.35 346
106.3 317
106.25 293
106.2 287
106.15 295
106.1 306
106.05 291
106 306
105.95 312
105.9 302
105.85 272
105.8 309
105.75 284
105.7 270
105.65 286
105.6 266
105.55 281
105.5 256
105.45 272
105.4 288
105.35 295
105.3 299
105.25 292
105.2 311
105.15 307
105.1 305
105.05 291
105 324
104.95 333
104.9 350
104.85 383
104.8 424
104.75 430
104.7 469
104.65 526
104.6 548
104.55 608
104.5 638
104.45 700
104.4 793
104.35 887
104.3 983
104.25 1041
104.2 1146
104.15 1222
104.1 1312
104.05 1389
104 1558
103.95 1602
103.9 1700
103.85 1840
103.8 1983
103.75 2073
103.7 2147
103.65 2257
103.6 2347
103.55 2370
103.5 2432
103.45 2484
103.4 2487
103.35 2473
103.3 2508
103.25 2497
103.2 2508
103.15 2454
103.1 2327
103.05 2289
103 2326
102.95 2241
102.9 2127
102.85 1871
102.8 1855
102.75 1685
102.7 1602
102.65 1453
102.6 1411
102.55 1300
102.5 1154
102.45 1067
102.4 987
102.35 934
102.3 842
102.25 790
102.2 742
102.15 683
102.1 644
102.05 629
102 573
101.95 569
101.9 579
101.85 567
101.8 536
101.75 524
101.7 521
101.65 508
101.6 517
101.55 514
101.5 478
101.45 440
101.4 440
101.35 485
101.3 467
101.25 460
101.2 468
101.15 483
101.1 485
101.05 509
101 494
100.95 480
100.9 490
100.85 487
100.8 485
100.75 473
100.7 476
100.65 456
100.6 436
100.55 462
100.5 506
100.45 559
100.4 661
100.35 729
100.3 872
100.25 998
100.2 1199
100.15 1579
100.1 1870
100.05 2273
100 2656
99.95 3019
99.9 3227
99.85 3339
99.8 3207
99.75 3017
99.7 2852
99.65 2683
99.6 2671
99.55 2893
99.5 3389
99.45 3916
99.4 4621
99.35 5264
99.3 5484
99.25 5662
99.2 5299
99.15 4596
99.1 3847
99.05 2929
99 2234
98.95 1572
98.9 1045
98.85 670
98.8 485
98.75 357
98.7 299
98.65 226
98.6 170
98.55 174
98.5 187
98.45 149
98.4 117
98.35 131
98.3 142
98.25 147
98.2 154
98.15 130
98.1 120
98.05 125
98 116
97.95 128
97.9 130
97.85 122
97.8 130
97.75 110
97.7 121
97.65 125
97.6 125
97.55 135
97.5 120
97.45 127
97.4 127
97.35 120
97.3 114
97.25 115
97.2 100
97.15 106
97.1 122
97.05 117
97 119
96.95 115
96.9 135
96.85 113

jeanbiego

2020/10/20 04:14

このままだとコードが読みづらいので、<code>ボタンを使ってインデントがわかるように編集してください。 https://teratail.com/help

y_waiwai

2020/10/20 04:15

このままではコードが読めないので、質問を編集し、<code>ボタンを押し、出てくる’’’の枠の中にコードを貼り付けてください

Yu_a

2020/10/20 04:24

お二方ともありがとうございます。

jeanbiego

2020/10/20 04:51

複雑なので私が回答できるかがまず不明ですが、それでも疑問がいくつか… 以下を質問文で明らかにしておいたほうが、有用な回答がつくのではないかと。・それぞれの変数定義は文献のどの式に対応していますか。・関数ではなくクラスで設計する意図はなんですか。・self.xの定義はどこですか。・yieldを使ってジェネレーターにしようとしているのは、あとでΣに入れるためですか。・「listの中に値が入らない」とありますが、動かすとどうなるのですか。・コピペで動かせるようにSi2p.txtか、それに代わるテストデータを貼れますか。

Yu_a

2020/10/21 04:37 編集

返信が遅くなり申し訳ありませんでした。その理由としては頂いた疑問点からプログラムを修正してみることに時間がかかってしまいました。自分がしたいことについてはPDFにまとめてみたので、それぞれの変数定義についてはわかると思います。listの中に値を入れるという点に対しても説明しておきました。Si2pのデータについても同様に示しておきます。どうぞよろしくお願いいたします。

jeanbiego

2020/10/21 05:07

まとめPDFはどこでしょうか。ちなみに、可能ならば質問文をそのまま編集してもらったほうが良いのですが。

jeanbiego

2020/10/21 05:14

失礼、すれ違いでしたね。 Q1 = h * (y1 - 0.5*(y[k] - y[-1] + y[203] - y[-1])) とかでy[-1]まあ最後の要素ということでわかります。y[203]は、yの203番めの要素ならy[202]としないといけませんが、これは要はy[-1]と同じです。この辺の意図はなんでしょうか。

jeanbiego

2020/10/21 05:19

Qの右辺の中の、(jx-b)+(jk-b)において、b=y[-1], jx=y[202], jk=y[k]で、この場合はb=jxってことですか？

Yu_a

2020/10/21 05:33

y[-1]とy[202]で分けてしまっているのはこちらのミスです。すみません。ｂ＝y[-1]=y[202], jx=y[k]です。ｋは図中ではデータの最終点ですが、プログラム内では任意の変数として使ってしまっていて混乱させてしまっています。

行動規範の内容に同意します

回答1件

ベストアンサー

while True:が無限ループなのでコメントアウト
y[203]をy[202]に修正
k, tはfor文でカウントアップ

とやるとグラフが出てきました。まだちゃんとPDFの中身を理解してないので、これが本当にご希望のものかは不明ですが、一応のせます。

python3
1import numpy as np
2import  matplotlib.pyplot as plt
3
4x,y = np.loadtxt ("Si2p.txt",skiprows=1,unpack = True)
5#shirley法によるバックグラウンドのプログラム
6Blist = []
7y1 = np.sum(y) - y[-1]
8h = 0.05
9
10
11Flist =[]
12#while True:
13for k in range(203):
14    R1 = h * (y1 - 0.5*(y[1] - y[-1]) + (y[k] - y[-1]))
15    Q1 = h * (y1 - 0.5*(y[k] - y[-1] + y[202] - y[-1]))
16    B1 = (y[1] - y[202])*(Q1/R1) + y[-1]
17    F = y[k] - B1
18    Blist.append(B1)
19    Flist.append(F)
20for t in range(100):
21    z = Flist[t]
22    Q = h * (y1 - 0.5*(z - y[-1] + y[202] - y[-1]))
23    R = h * (y1 - 0.5*(y[1] - y[-1]) + (z - y[-1]))
24    B = (y[1] - y[202])*(Q/R) + y[-1]
25    Blist.append(B)
26    Flist.append(F)
27    
28#print(Blist)
29
30plt.scatter(x,y,marker="o" ,color="black",s=10)
31ax = plt.gca()        
32y_min, y_max = ax.get_ylim()
33ax.set_ylim(0,6000)
34x_min, x_max = ax.get_xlim()
35#ax.set_xlim(97,102)
36plt.title("Shirley",fontname="MS Gothic")
37plt.plot(x,y,color="black")
38plt.xlabel("BandEnergy[eV]",fontname="MS Gothic")
39plt.ylabel("Intensity[a.u.]",fontname="MS Gothic")
40ax.invert_xaxis()
41
42plt.show()
43plt.savefig("Shirley.png")