F(x=y/a)な関数を基に、F1(y=a1x),F2(y=a2x) ~みたいに変換したい

原型のF(x=s/a)から、a1,a2,a3 ~ の定数項で変化を持たせた F1(y) F2(y) F3(y) ~ を作りたいけど方法がわからない

原型の関数が下のコード

python
1from scipy.special import kv
2import matplotlib.pyplot as plt
3from scipy.integrate import quad
4import numpy as np
5import math
6from math import gamma
7
8xs = np.logspace(-4, 2, 10000)
9
10f = lambda z: kv(5/3,z)
11F = [quad(f,x,np.inf)[0]*x for x in xs]
12
13a = gamma(1/3) 
14G = [(4*math.pi/np.sqrt(3)/a)*(x/2)**(1/3) for x in xs]
15H = [((math.pi/2)**(1/2))*(x**(1/2))*(math.exp(-x))for x in xs]
16
17
18
19def A(x,F,G,H):
20    if x <= 5.0*1e-3:
21        g =G
22    elif 5.0*1e-3 < x < 30:
23        g =F
24    elif 30 <= x:
25        g =H
26    return g
27
28vec_A = np.vectorize(A)
29y = vec_A(xs,F,G,H)
30

ここから、x→ s/a1,s/a2 ~ みたいに変化させたいけど書き方の見当がつかない

python
1from scipy.special import kv
2import matplotlib.pyplot as plt
3from scipy.integrate import quad
4import numpy as np
5import math
6from math import gamma
7
8m_e=9.10938356*1e-31
9c=2.99792458*1e+8
10q=1.6021766208*1e-19
11
12pa=math.pi/2
13sin_pa=math.sin(pa)
14from sympy import *
15
16lf1=10
17lf2=100
18lf3=1000
19
20B1=10*1e-6*1e-4
21B2=100*1e-6*1e-4
22
23omega_c1=3*lf1**2*q*B1*sin_pa/(2*m_e*c)
24omega_c2=3*lf2**2*q*B1*sin_pa/(2*m_e*c)
25omega_c3=3*lf3**2*q*B1*sin_pa/(2*m_e*c)
26omega_c4=3*lf1**2*q*B2*sin_pa/(2*m_e*c)
27omega_c5=3*lf2**2*q*B2*sin_pa/(2*m_e*c)
28omega_c6=3*lf3**2*q*B2*sin_pa/(2*m_e*c)
29
30###################################################1#
31###x→ω x=x*ω_c#
32xs1 = np.linspace(1e-4*omega_c1, 1e2*omega_c1, 10000)
33
34
35f = lambda z: kv(5/3,z)
36F = [quad(f,x/omega_c1,np.inf)[0]*(x/omega_c1) for x in xs1]
37
38a = gamma(1/3) 
39G = [(4*math.pi/np.sqrt(3)/a)*((x/omega_c1)/2)**(1/3) for x in xs1]
40H = [((math.pi/2)**(1/2))*((x/omega_c1)**(1/2))*(math.exp(-(x/omega_c1)))for x in xs1]
41
42
43
44def A(x,F,G,H):
45    if (x/omega_c1) <= 5.0*1e-3:
46        g =G
47    elif 5.0*1e-3 < (x/omega_c1) < 30:
48        g =F
49    elif 30 <= (x/omega_c1):
50        g =H
51    return g
52
53vec_A = np.vectorize(A)
54y = vec_A(xs1,F,G,H)
55
56P1=(math.sqrt(3)*q**3*B1*sin_pa/(2*math.pi*m_e*c**2))*y
57###################################################2#
58xs2 = np.linspace(1e-4*omega_c2, 1e2*omega_c2, 10000)
59
60
61f = lambda z: kv(5/3,z)
62F = [quad(f,x/omega_c2,np.inf)[0]*(x/omega_c2) for x in xs2]
63
64a = gamma(1/3) 
65G = [(4*math.pi/np.sqrt(3)/a)*((x/omega_c2)/2)**(1/3) for x in xs2]
66H = [((math.pi/2)**(1/2))*((x/omega_c2)**(1/2))*(math.exp(-(x/omega_c2)))for x in xs2]
67
68
69
70def A(x,F,G,H):
71    if (x/omega_c2) <= 5.0*1e-3:
72        g =G
73    elif 5.0*1e-3 < (x/omega_c2) < 30:
74        g =F
75    elif 30 <= (x/omega_c2):
76        g =H
77    return g
78
79vec_A = np.vectorize(A)
80y = vec_A(xs2,F,G,H)
81
82P2=(math.sqrt(3)*q**3*B1*sin_pa/(2*math.pi*m_e*c**2))*y
83###################################################3#
84xs3 = np.linspace(1e-4*omega_c3, 1e2*omega_c3, 10000)
85
86
87f = lambda z: kv(5/3,z)
88F = [quad(f,x/omega_c3,np.inf)[0]*(x/omega_c3) for x in xs3]
89
90a = gamma(1/3) 
91G = [(4*math.pi/np.sqrt(3)/a)*((x/omega_c3)/2)**(1/3) for x in xs3]
92H = [((math.pi/2)**(1/2))*((x/omega_c3)**(1/2))*(math.exp(-(x/omega_c3)))for x in xs3]
93
94
95
96def A(x,F,G,H):
97    if (x/omega_c3) <= 5.0*1e-3:
98        g =G
99    elif 5.0*1e-3 < (x/omega_c3) < 30:
100        g =F
101    elif 30 <= (x/omega_c3):
102        g =H
103    return g
104
105vec_A = np.vectorize(A)
106y = vec_A(xs3,F,G,H)
107
108P3=(math.sqrt(3)*q**3*B1*sin_pa/(2*math.pi*m_e*c**2))*y
109###################################################4#
110xs4 = np.linspace(1e-4*omega_c4, 1e2*omega_c4, 10000)
111
112
113f = lambda z: kv(5/3,z)
114F = [quad(f,x/omega_c4,np.inf)[0]*(x/omega_c4) for x in xs4]
115
116a = gamma(1/3) 
117G = [(4*math.pi/np.sqrt(3)/a)*((x/omega_c4)/2)**(1/3) for x in xs4]
118H = [((math.pi/2)**(1/2))*((x/omega_c4)**(1/2))*(math.exp(-(x/omega_c4)))for x in xs4]
119
120
121
122def A(x,F,G,H):
123    if (x/omega_c4) <= 5.0*1e-3:
124        g =G
125    elif 5.0*1e-3 < (x/omega_c4) < 30:
126        g =F
127    elif 30 <= (x/omega_c4):
128        g =H
129    return g
130
131vec_A = np.vectorize(A)
132y = vec_A(xs4,F,G,H)
133
134P4=(math.sqrt(3)*q**3*B2*sin_pa/(2*math.pi*m_e*c**2))*y
135###################################################5#
136xs5 = np.linspace(1e-4*omega_c5, 1e2*omega_c5, 10000)
137
138
139f = lambda z: kv(5/3,z)
140F = [quad(f,x/omega_c5,np.inf)[0]*(x/omega_c5) for x in xs5]
141
142a = gamma(1/3) 
143G = [(4*math.pi/np.sqrt(3)/a)*((x/omega_c5)/2)**(1/3) for x in xs5]
144H = [((math.pi/2)**(1/2))*((x/omega_c5)**(1/2))*(math.exp(-(x/omega_c5)))for x in xs5]
145
146
147
148def A(x,F,G,H):
149    if (x/omega_c5) <= 5.0*1e-3:
150        g =G
151    elif 5.0*1e-3 < (x/omega_c5) < 30:
152        g =F
153    elif 30 <= (x/omega_c5):
154        g =H
155    return g
156
157vec_A = np.vectorize(A)
158y = vec_A(xs5,F,G,H)
159
160P5=(math.sqrt(3)*q**3*B2*sin_pa/(2*math.pi*m_e*c**2))*y
161###################################################6#
162xs6= np.linspace(1e-4*omega_c6, 1e2*omega_c6, 10000)
163
164
165f = lambda z: kv(5/3,z)
166F = [quad(f,x/omega_c6,np.inf)[0]*(x/omega_c6) for x in xs6]
167
168a = gamma(1/3) 
169G = [(4*math.pi/np.sqrt(3)/a)*((x/omega_c6)/2)**(1/3) for x in xs6]
170H = [((math.pi/2)**(1/2))*((x/omega_c6)**(1/2))*(math.exp(-(x/omega_c6)))for x in xs6]
171
172
173
174def A(x,F,G,H):
175    if (x/omega_c6) <= 5.0*1e-3:
176        g =G
177    elif 5.0*1e-3 < (x/omega_c6) < 30:
178        g =F
179    elif 30 <= (x/omega_c6):
180        g =H
181    return g
182
183vec_A = np.vectorize(A)
184y = vec_A(xs6,F,G,H)
185
186P6=(math.sqrt(3)*q**3*B2*sin_pa/(2*math.pi*m_e*c**2))*y
187
188#####################################################
189
190fig = plt.figure()
191ax = fig.add_subplot(1,1,1)
192ax.grid()
193ax.set_xlim(1e-4*omega_c1,  10000*1e2*omega_c1)
194ax.set_ylim(1e-3*(math.sqrt(3)*q**3*B1*sin_pa/(2*math.pi*m_e*c**2)), 100*1e0*(math.sqrt(3)*q**3*B1*sin_pa/(2*math.pi*m_e*c**2)))
195ax.set_yscale('log')  
196ax.set_xscale('log')  
197ax.set_title('')
198ax.set_xlabel('ω')
199ax.set_ylabel('P')
200ax.plot(xs1,P1,"c",label="1")
201ax.plot(xs2,P2,"r",label="2")
202ax.plot(xs3,P3,"g",label="3")
203ax.plot(xs4,P4,"b",label="4")
204ax.plot(xs5,P5,"black",label="5")
205ax.plot(xs6,P6,"y",label="6")
206ax.legend()
207
208plt.show