1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91
| #!/usr/bin/env python
import math
import matplotlib.pyplot as plt
import numpy as np
import sys
# Unites CGS
G = 6.6726e-8
sigma = 5.6705e-5
PI = math.pi
sys.path.append('/home/damien/COURS/PIR/ester/python/')
from star import *
indice = 0
puiss = bucket = [0] * 1500
masse = bucket = [0] * 1500
omega = bucket = [0] * 1500
compo = bucket = [0] * 1500
# Boucle sur les fichiers
# print 'je suis devant la boucle'
for i in range(30,35,5) :
for j in range(2,4,1) :
for k in range(2,4,1) :
# print 'je suis rentre dans la boucle'
if i > 50 :
k=k*10
# print i
# print j
# print k
mass = str(i)
vitess = str(j)
comp = str(k)
seq = [mass,vitess,comp]
A = star2d("-".join(seq))
Re = A.Re
M = A.M
L = A.L
theta = A.th
# Constantes de normalisation de Teff et geff
C1 = pow(L/(4*PI*sigma*Re*Re),-0.25)
C2 = math.pow(G*M/(Re*Re),-1)
# Teff et geff + en log pour l'affichage
teff = C1*(A.Teff)
geff = C2*(A.gsup)
tefflog = np.log10(teff)
gefflog = np.log10(geff)
# Fit a=coeff directeur de la droite = puissance entre Teff et geff
a,b = np.polyfit(gefflog, tefflog, 1)
y_new = a*gefflog +b
c,d,e = np.polyfit(gefflog, tefflog, 2)
y_new2 = c*gefflog*gefflog + d*gefflog + e
puiss[indice] = a
if c < a/10 :
puiss[indice] = d
masse[indice] = i/10.0
omega[indice] = j/10.0
compo[indice] = k/10.0
if i > 50 :
compo[indice] = k/100
indice = indice + 1
# Affichage
#plt.plot(np.array(gefflog),np.array(tefflog))
#plt.plot(gefflog,y_new)
#plt.plot(gefflog,y_new2)
#plt.xlabel('log10(geff) normalisee')
#plt.ylabel('log10(Teff) normalisee')
#plt.show() |
Partager