int r = 3;
int v = 5;
int x = 400;
int y = 400;
int c=0;
void setup() {
size(1000, 1000);
background(255);
noStroke();
fill(c);
x=width/2;
y=height/2;
}
void draw() {
fill(c);
float th=random(0,2*PI);
x += int(v*cos(th));
y += int(v*sin(th));
if(x>width) x=0;
if(x<0 if="" x="width;" y="">height) y=0;
if(y<0 c="" if="" y="height;">255) c=0;
}
import processing.serial.*;
Serial port;
PImage img;
float Y=0;
void setup() {
size(640, 640, P3D);
img=loadImage("Kao.jpg");
port = new Serial(this,"COM3",9600);
frameRate(30);
stroke(0, 0, 0, 100);
//port.bufferUntil('\n');
port.write('a');
}
void draw() {
lights();
background(0);
float cameraY = height/2.0;
float fov = mouseX/float(width) * PI/6;
float cameraZ = cameraY / tan(fov / 2.0);
float aspect = float(width)/float(height);
perspective(fov, aspect, cameraZ/10.0, cameraZ*10.0);
translate(width/2, height/2, 0);
rotateX((1.9- Y/9.6) * PI/2);
rotateZ(PI/3 + mouseX/float(height) * PI);
beginShape();
box(40,40,50);
translate(10, 0, -50);
box(15,15,40);
translate(-20, 0, 0);
box(15,15,40);
translate(-20, 0, 50);
box(15,15,40);
translate(60, 0, 0);
box(15,15,40);
translate(-30, 0, 40);
box(25,25,25);
translate(0, -12.5, 0);
box(25,1,25);
texture(img);
vertex(0,100,0,img.width,img.height);
endShape();
translate(0, -10, 20);
rotateX(radians(-90));
rotateY(radians(180));
if(Y<4)
{
text("H!",0,0,0);
}
else if(mousePressed)
{
text("I Love You!",0,0,0);
}
if(port.available()>1){
String ms=port.readStringUntil('\n');
if(ms!=null){
ms=trim(ms);
float Y0=float(ms);
Y=int(Y0*10)/10;
}
port.write("a");
}
print(Y);
print("\n");
}
void mousePressed(){
port.write('a');
}
Arudino側。
こちら
https://github.com/ArduSat/ArdusatSDK
を利用。
#include Arduino.h
#include Wire.h
#include ArdusatSDK.h
ArdusatSerial serialConnection(SERIAL_MODE_HARDWARE_AND_SOFTWARE, 8, 9);
Acceleration accel;
void setup(void)
{
Serial.begin(9600);
serialConnection.println("a");
}
void loop(void)
{
int a,c;
if(Serial.available()>0){
accel.readToJSON("accelerometer");
c=Serial.read();
Serial.println(accel.z);
}
}
PImage img;
void setup() {
size(640, 640, P3D);
img=loadImage("Kao.jpg");
stroke(0, 0, 0, 100);
}
void draw() {
lights();
background(0);
float cameraY = height/2.0;
float fov = mouseX/float(width) * PI/6;
float cameraZ = cameraY / tan(fov / 2.0);
float aspect = float(width)/float(height);
perspective(fov, aspect, cameraZ/10.0, cameraZ*10.0);
translate(width/2, height/2, 0);
rotateX(PI/3 + mouseY/float(height) * PI);
rotateZ(PI/3 + mouseX/float(height) * PI);
beginShape();
box(40,40,50);
translate(10, 0, -50);
box(15,15,40);
translate(-20, 0, 0);
box(15,15,40);
translate(-20, 0, 50);
box(15,15,40);
translate(60, 0, 0);
box(15,15,40);
translate(-30, 0, 40);
box(25,25,25);
translate(0, -12.5, 0);
box(25,1,25);
texture(img);
vertex(0,100,0,img.width,img.height);
endShape();
translate(0, -10, 20);
rotateX(radians(-90));
rotateY(radians(180));
if(PI/3 + mouseY/float(height) * PI>2.3){
text("H!",0,0,0);
}
else if(mousePressed) {
text("Love :-)",0,0,0);
}
}
# Japanese is written in UTF-8.
# シャープはコメントです。
##勉強用制御
No1=0
No2=0
No3=0
##色々設定
bin=50 # bin面体サイコロ
dice_test=10000 # サイコロのチェックのために振る回数
one_test=10 # one_test回振って平均を求める
test=1000 # 平均を求める試行回数
##出力ファイル
dice_check="dice_check_bin${bin}_test${diec_test}"
dice_check_txt="${dice_check}.txt"
dice_check_ps="${dice_check}.ps"
CLT_check="CLT_check_bin${bin}_test${test}"
CLT_check_txt="${CLT_check}.txt"
CLT_check_ps="${CLT_check}.ps"
rm -f $dice_check_txt
rm -f $dice_check_ps
rm -f $CLT_check_txt
rm -f $CLT_check_ps
#さいころ関数
dice_func(){
echo "$RANDOM" $bin | awk '{print int($1/32767*$2)}'
}
#サイコロ関数のテスト
echo "Dice is" `dice_func 0`"."
####
if (($No1)) ;
then
####
#配列の初期化1
array=()
for((i=0;i do
array[${i}]=0
done
#サイコロ関数の正しさをtest回ふって確認。
for((i=0;i do
rand=`dice_func 0`
array[${rand}]=$((${array[${rand}]}+1))
done
#結果をファイルに書き込み
for((i=0;i do
echo "$i ${array[${i}]}" >> ${dice_check_txt}
done
echo ${dice_check_txt} "is made."
###
fi #No1
###
####
if (($No2)) ;
then
####
#配列の初期化2
for((i=0;i do
array[${i}]=0
done
#「サイコロ関数をone_test回ふってその平均を調べる。」試行をtest回実施。
for((i=0;i do
s=0
for((j=0;j do
rand=`dice_func 0`
s=$((${s}+${rand}))
done
ave=`echo $s $one_test | awk '{print int($1/$2)}'`
array[${ave}]=$((${array[${ave}]}+1))
echo $i "kai end."
done
#結果をファイルに書き込み
for((i=0;i do
echo "$i ${array[${i}]}" >> ${CLT_check_txt}
done
echo ${CLT_check_txt} "is made."
###
fi #No2
###
####
if (($No3)) ;
then
####
#初期値を探す
n0=`awk 'BEGIN{max_n=0}{if($2>max_n) max_n=$2}END{print max_n}' ${CLT_check_txt}`
m0=`awk 'BEGIN{max_m=0;max_n=0}{if($2>max_n){max_n=$2;max_m=$1}}END{print max_m}' ${CLT_check_txt}`
#Gnuplotで作図
gnuplot<
set ylabel "Number"
set xlabel "Value"
set te po co solid
set title "Dice Check : Bin ${bin} Test ${dice_test}"
set ou "${dice_check_ps}"
pl [][0:] "${dice_check_txt}" w histeps lw 2
set ou
Gauss(x)=n*exp(-0.5*((x-m)/s)**2)
n=${n0}
m=${m0}
s=m*0.1
fit Gauss(x) "${CLT_check_txt}" via n,m,s
set title "CLT Check : Ave. of ${one_test} Bin ${bin} (Test ${dice_test})"
set ou "${CLT_check_ps}"
pl "${CLT_check_txt}" w histeps lw 2, Gauss(x)
set ou
exit
EOF
###
fi #No3
###
T=int(1e18) #これが1e9だと不可。うーん...。
def cut0(val) :
while True :
if val%10==0 :
val=val/10
else :
break
return val
def func1(N) :
val=1
for num in xrange(1,N+1) :
val=cut0(num*val)
val=val%T
return val%int(1e9)
print func1(int(raw_input()))
#! /bin/bash echo "Inu" echo "Saru" echo "Kiji" :<<"_Test_" echo "Neko" ls xspec<<EOF exit EOF _TEST_
SCRIPT_SHELL="python"
SCRIPT_NAME="Arf_BG_corrector.py"
SCRIPT_HELP_ARG="1) Input Rate Bg PI file made by mathpha 2) Src Arf 3) Bg Arf 4) Output Bg PI file"
SCRIPT_HELP_SENTENCE="Modiefy Bg PI Channel with Src / Bg arf ratio."
SCRIPT_NUM_ARG=4
###IMPORT###
import sys
import pyfits
import os
###HELP###
if len(sys.argv[1:])!=SCRIPT_NUM_ARG:
print 'Name: '+SCRIPT_NAME
print 'Arguments: '+SCRIPT_HELP_ARG
print 'Explanation: '+SCRIPT_HELP_SENTENCE
sys.exit()
###MAIN###
ev_per_ch=3.65
bg_pif=str(sys.argv[1])
src_arf=str(sys.argv[2])
bg_arf=str(sys.argv[3])
out_pif=str(sys.argv[4])
def main() :
if os.path.isfile(out_pif) :
os.remove(out_pif)
data_src_arf = (pyfits.open(src_arf))['SPECRESP'].data
data_bg_arf = (pyfits.open(bg_arf))['SPECRESP'].data
ratio=[]
for n in range(len(data_src_arf)) :
ratio.append([data_src_arf['ENERG_LO'][n],data_src_arf['ENERG_HI'][n],data_src_arf['SPECRESP'][n] / data_bg_arf['SPECRESP'][n]])
hdulist=pyfits.open(bg_pif)
sp=hdulist['SPECTRUM'].data
ch_num=len(hdulist['SPECTRUM'].data)
hdulist['SPECTRUM'].header['COMMENT']='Modefied with Arf_BG_corrector.py'
hdulist['SPECTRUM'].header['COMMENT']='Src arf : '+src_arf
hdulist['SPECTRUM'].header['COMMENT']='Bg arf : '+bg_arf
hdulist['SPECTRUM'].header['COMMENT']='Bg pi : '+bg_pif
for n in xrange(ch_num) :
enegy=ev_per_ch*n/1e3
for erange in ratio :
if erange[0] <= enegy and enegy < erange[1] :
sp['RATE'][n]=sp['RATE'][n]*erange[2]
sp['STAT_ERR'][n]=sp['STAT_ERR'][n]*erange[2]
break
hdulist.writeto(out_pif)
main()
###IMPORT###
import sys
import openpyxl
from openpyxl.cell import get_column_letter
###MAIN###
def main() :
inbook_file="Book1.xlsx"
csv_file="CSV.csv"
out_file="Book2.xlsx"
#Read XLSX
wb1=openpyxl.load_workbook(inbook_file)
ws1=wb1.active
AB={}
for u in ws1.rows :
AB[str(u[2].value)]=[u[0].value,u[1].value]
#Read CSV
AL=[]
for line in open(csv_file, 'r'):
AL.append((line.split(","))[0])
#Open WB for Output
wb2=openpyxl.workbook.Workbook()
ws2=wb2.active
for i,ad in enumerate(AL) :
ws2['A'+str(i+1)].value=ad
if AB.has_key(ad) :
ws2['B'+str(i+1)]=AB[ad][0]
ws2['C'+str(i+1)]=AB[ad][1]
else :
ws2['B'+str(i+1)]=" "
ws2['C'+str(i+1)]=" "
wb2.save(filename=out_file)
main()
SCRIPT_SHELL="python"
SCRIPT_NAME="PiShift.py"
SCRIPT_HELP_ARG="1) input Count PI file made by mathpha 2)output PI file 3)c0 4)c1 5)c2"
SCRIPT_HELP_SENTENCE="Modify PI Channel with f(x)=c2*x**2+c1*x+c0. Only for Suzaku XIS."
SCRIPT_NUM_ARG=5
###IMPORT###
import sys
import pyfits
import os
###HELP###
if len(sys.argv[1:])!=SCRIPT_NUM_ARG:
print 'Name: '+SCRIPT_NAME
print 'Arguments: '+SCRIPT_HELP_ARG
print 'Explanation: '+SCRIPT_HELP_SENTENCE
sys.exit()
###MAIN###
inf=str(sys.argv[1])
outf=str(sys.argv[2])
c0=float(sys.argv[3])
c1=float(sys.argv[4])
c2=float(sys.argv[5])
ch_num=4096
hdulist = pyfits.open(inf)
if os.path.isfile(outf) :
os.remove(outf)
sp=hdulist[1].data
hdulist[1].header['COMMENT']='Modefied with PiSift.py'
hdulist[1].header['COMMENT']='c0='+str(c0)
hdulist[1].header['COMMENT']='c1='+str(c1)
hdulist[1].header['COMMENT']='c2='+str(c2)
def f(x) :
x=float(x)
return c2*x**2+c1*x+c0
def main() :
new_counts=[]
for n in range(ch_num) :
new_counts.append(0.0)
old_counts=sp['COUNTS']
for n in range(ch_num) :
new_ch0=f(n)
new_ch1=f(n+1)
if new_ch1 < new_ch0 :
sys.exit('Correction function is not monotone.')
elif new_ch0 > 0 and new_ch1 < ch_num :
for m in range(int(new_ch0+1),int(new_ch1)) :
new_counts[m]+=old_counts[n]/(new_ch1-new_ch0)
new_counts[int(new_ch0)]+=old_counts[n]/(new_ch1-new_ch0)*(int(new_ch0+1)-new_ch0)
new_counts[int(new_ch1)]+=old_counts[n]/(new_ch1-new_ch0)*(new_ch1-int(new_ch1))
for n in range(ch_num) :
sp['COUNTS'][n]=int(round(new_counts[n]))
sp['STAT_ERR'][n]=float(new_counts[n]**0.5)
hdulist.writeto(outf)
main()
data=map(lambda str:str.split(), open("FeDist_Screened.plt").readlines())
とはいえ、これはメンドくさすぎると感じるのは自分だけ?dtypeどうにかならんのか。 dtype=('f8','f8','S16')とかダメなのか...。
data=np.loadtxt('FeDist_Screened.plt', delimiter=' ', usecols=[0,2,11],skiprows=1,dtype={'names':('l','6.4','id'), 'formats':('f8','f8','S16')})
一応やりたかったことをやってみた結果。awkだと瞬間なんだけど...。
SCRIPT_TYPE="InstrumentalityOfMankind"
SCRIPT_SHELL="python"
SCRIPT_NAME="make_Data.py"
SCRIPT_HELP_ARG=""
SCRIPT_HELP_SENTENCE=""
SCRIPT_NUM_ARG=0
SCRIPT_VERSION=1.0
###IMPORT###
import sys
###HELP###
if len(sys.argv[1:])!=SCRIPT_NUM_ARG:
print 'Name: '+SCRIPT_NAME
print 'Arguments: '+SCRIPT_HELP_ARG
print 'Explanation: '+SCRIPT_HELP_SENTENCE
sys.exit()
###MAIN###
infile="FeDist_Screened.plt"
data=map(lambda str:str.split(), open("FeDist_Screened.plt","r").readlines())
f=1/(3.14*18**2)
enes=['6.4','6.7']
for ene in enes :
if ene=="6.4" :
line_num=2
elif ene=="6.7" :
line_num=5
outname='Fe'+ene+'Dist_GRXE.txt'
line_emin_num=line_num+1
line_emax_num=line_num+2
outfile=open(outname,"w")
for data_row in data[1:] :
if float(data_row[0])**2>10**2:
l,b=data_row[0:2]
ids=data_row[11]+' '+data_row[12]
val, val_min, val_max=map(float,data_row[line_num:line_num+3])
val=f*val
val_err=f*(val_max-val_min)*0.5/1.64
out_str=l+' '+b+' '+str(val)+' '+str(val_err)+' '+ids+'\n'
outfile.write(out_str)
outfile.close()
書き込む前にstrにするのとかが果てしなく面倒くさい。
結論 : awkで良いじゃん (ヲイ) or 適材適所。
###IMPORT###
import sys
import random
###MAIN###
def dice():
return int(random.uniform(0,365))
def main():
trial=1e5
count=[]
same=0
N=40
M=3
for day in range(0,365) :
count.append(0)
for num1 in range(0,int(trial)):
for day in range(0,365) :
count[day]=0
for num2 in range(0,N) :
count[dice()]+=1
for day in range(0,365) :
if count[day] > (M-1) :
same+=1
break
print same/trial
main()
###MAIN###
N=5
def perceptron(m,x,w,v):
prec=[]
if m<=0:
return -1
for k in range(0,3) :
prec.append(-1)
prec[m-1]=1
for k in range(0,3):
if prec[k]*v[k]<=0 :
for j in range(0,3*N+1) :
w[(3*N+1)*k+j]+=prec[k]*x[j]
for i in range(0,3*N-3) :
x[3*N-1-i]=x[3*N-4-i]
for i in range(0,3):
x[i]=prec[i]
for k in range(0,3):
v[k]=0
for k in range(0,3):
for j in range(0,3*N+1):
v[k]+=w[(3*N+1)*k+j]*x[j]
vmax=-1000000
for k in range(0,3) :
if v[k] >=vmax :
vmax=v[k]
kmax=k
return kmax+1
def main() :
total=0
v=[]
fw=[]
x=[]
w=[]
pred=0
for i in range(0,3):
v.append(0)
fw.append(0)
for i in range(0,3*N):
x.append(0)
x.append(-1)
for i in range(0, 9*N+3):
w.append(0)
print "x:",x
print "w:",w
print "v:",v
m=1
while m>0:
pred=perceptron(m,x,w,v)
while True:
try:
m = int(raw_input("1 Gu 2 Choki 3 Pa :"))
break
except ValueError:
print "Try again..."
if m>3 :
m=3
print "You",m,"vs Machine",(pred+1)%3+1
if pred==m :
print "Machine win."
fw[0]+=1
elif pred%3==m-1:
print "You win"
fw[1]+=1
else:
print "Draw"
fw[2]+=1
total+=1
print "You :",fw[1], "Machine:", fw[0], "Draw :",fw[2],"Total:", total
#print "x:",x
#print "w:",w
#print "v:",v
return 0
main()
import random
def main():
vals=['ST','CON','LK','DEX','IQ','CHR','SPD','MON']
ch={}
for val in vals:
ch[val]=int(three_dices())
bonus=cal_bonus(ch)
print ch,bonus
def cal_bonus(ch) :
vals=['ST','LK','DEX']
bonus=0
for val in vals :
if ch[val]>12 :
bonus=bonus+ch[val]-12
elif ch[val]<9 :
bonus=bonus+ch[val]-9
return bonus
def dice():
random.seed()
return random.randint(1, 6)
def three_dices():
s=0
for num in [0,1,2]:
s=s+dice()
return s
main()
import random
import numpy
import matplotlib.pyplot as plt
def normal_dice(dice_num) :
vals = []
total_val = 0
random.seed()
for num in range(0,dice_num):
vals.append(random.randint(1, 6))
#print vals
total_val = sum(vals)
return total_val
def berserk_dice(dice_num) :
vals = []
total_val = 0
random.seed()
for num in range(0,dice_num):
vals.append(random.randint(1, 6))
#print vals
total_val = total_val + sum(vals)
for num in range(1,6):
num_of_same = vals.count(num)
if num_of_same > 1 :
#print "Bersrk by %d dices of %d marks" % (num_of_same,num)
total_val=total_val+berserk_dice(num_of_same)
return total_val
def hyper_berserk_dice(dice_num) :
vals = []
total_val = 0
random.seed()
for num in range(0,dice_num):
vals.append(random.randint(1, 6))
#print vals
total_val = total_val + sum(vals)
for num in range(1,6):
num_of_same = vals.count(num)
if num_of_same > 1 :
# print "Bersrk by %d dices of %d marks" % (num_of_same,num)
total_val=total_val+hyper_berserk_dice(num_of_same+1)
return total_val
def main():
dice_num = input("Input dice_num :")
trial = input("Input Num. of Trials :")
vals=[]
for num in range(0, trial) :
vals.append(berserk_dice(dice_num))
if num%10000==0:
print num
vals=numpy.array(vals)
print "Mean : ", numpy.mean(vals)
print "SD : ", numpy.std(vals)
print "Max : ", vals.max()
print "Min : ", vals.min()
plt.rc('font', **{'family': 'serif'})
fig = plt.figure()
ax = fig.add_subplot(111)
ax.hist(vals, bins=100)
ax.set_title('Berserk with %d dices! (%d trials)' % (dice_num, trial))
ax.set_xlabel('Hits!')
ax.set_ylabel('Frequency')
plt.show()
return None
main()
# -*- coding: utf-8 -*-
import random
import numpy
import matplotlib.pyplot as plt
'''
あまり格好は良くないけど、中心極限定理をテストするスクリプト。
歪んだベータ分布から正規分布が出てくる。
以下を参考にした。
http://akiyoko.hatenablog.jp/entry/2013/06/07/213448
http://docs.python.jp/2/library/random.html
'''
def main():
random_val=[]
print("This is a test of central limit theorem by beta functions.")
alpha = input("Input alpha :")
beta = input("Input beta :")
N = input("Input N :")
trial = input("Input Num. of Trials :")
random.seed()
for num in range(0, N):
random_val.append(random.betavariate(alpha,beta))
random_val=numpy.array(random_val)
mean=numpy.mean(random_val)
print "mean=%.3f" % mean
sigma=numpy.std(random_val)
print "sigma=%.3f" % sigma
vals=[]
for num0 in range(0,trial):
val=0.0
for num1 in range(0, N):
val=val+random.betavariate(alpha,beta)
vals.append(val)
vals=numpy.array(vals)
plt.rc('font', **{'family': 'serif'})
fig = plt.figure()
ax0 = fig.add_subplot(311)
ax0.hist(random_val, bins=100, range=(0, 1), normed=False, facecolor='r', alpha=0.8)
ax0.set_xlim(0, 1)
y,x=numpy.histogram(random_val, bins=100, range=(0, 1), normed=False)
ax0.set_title('Beta distribution')
ax0.set_xlabel('Random variable')
ax0.set_ylabel('Frequency')
ax0.text(0.8, y.max()*0.8,r'alpha=%.2f' % alpha)
ax0.text(0.8, y.max()*0.7,r'beta=%.2f' % beta)
ax0.text(0.8, y.max()*0.6,r'N=%d' % N)
ax1 = fig.add_subplot(313)
ax1.hist(vals, bins=100, normed=False, facecolor='g', alpha=0.8)
ax1.set_title('Distribution of sums Trial=%d' % trial)
ax1.set_xlabel('Sum of random variables')
ax1.set_ylabel('Frequency')
plt.show()
return None
main()
#!/Library/Frameworks/Python.framework/Versions/7.2/bin/python
# -*- coding: utf-8 -*-
SCRIPT_TYPE="InstrumentalityOfMankind"
SCRIPT_SHELL="python"
SCRIPT_NAME=""
SCRIPT_HELP_ARG=""
SCRIPT_HELP_SENTENCE=""
SCRIPT_NUM_ARG=0
SCRIPT_VERSION=1.0
###IMPORT###
import sys
import tweepy
import random
import time
###HELP###
if len(sys.argv[1:])!=SCRIPT_NUM_ARG:
print 'Name: '+SCRIPT_NAME
print 'Arguments: '+SCRIPT_HELP_ARG
print 'Explanation: '+SCRIPT_HELP_SENTENCE
sys.exit()
###MAIN###
CONSUMER_KEY = ''
CONSUMER_SECRET = ''
logfile="temp.log"
###ふぁんくしょん###
def Tweets(ack, acs, text) :
CONSUMER_KEY = ''
CONSUMER_SECRET = ''
auth = tweepy.OAuthHandler(CONSUMER_KEY, CONSUMER_SECRET)
auth.set_access_token(ack, acs)
api = tweepy.API(auth, api_root='/1.1')
api.update_status(text)
for p in tweepy.Cursor(api.user_timeline).items(1):
lastid=p.id
print text
del api
return lastid
def Tweets_as_A(text) :
ACCESS_KEY1 = ''
ACCESS_SECRET1 = ''
return Tweets(ACCESS_KEY1, ACCESS_SECRET1, text)
def Tweets_as_B(text) :
ACCESS_KEY2 = ''
ACCESS_SECRET2 = ''
return Tweets(ACCESS_KEY2, ACCESS_SECRET2, text)
def Reply(ack, acs, text, lastid, to_usr) :
CONSUMER_KEY = ''
CONSUMER_SECRET = ''
auth = tweepy.OAuthHandler(CONSUMER_KEY, CONSUMER_SECRET)
auth.set_access_token(ack, acs)
api = tweepy.API(auth, api_root='/1.1')
name_to_reply = to_usr
status_id_to_reply = lastid
text2 = u'@%s %s' % (name_to_reply, text)
api.update_status(text2, in_reply_to_status_id = status_id_to_reply)
print text2
for p in tweepy.Cursor(api.user_timeline).items(1):
a=p.id
del api
return a
def Reply_as_A(text,lastid) :
ACCESS_KEY1 = ''
ACCESS_SECRET1 = ''
return Reply(ACCESS_KEY1, ACCESS_SECRET1, text, lastid, "B")
def Reply_as_B(text, lastid) :
ACCESS_KEY2 = ''
ACCESS_SECRET2 = ''
return Reply(ACCESS_KEY2, ACCESS_SECRET2, text, lastid, "A")
####おだいけってい###
num_wadai=2
f1=open(logfile,"r")
old_nums=[]
for line in f1:
old_nums.append(line)
f1.close()
prev_num=int(old_nums[-1])
num_now=random.randint(1, num_wadai)
while prev_num == num_now :
num_now=random.randint(1, num_wadai)
print num_now, prev_num
f1=open(logfile,"w")
old_nums.append(str(num_now)+"\n")
f1.writelines(old_nums)
f1.close()
####おしゃべり#####
if num_now==1 :
text=u"はーちみつ！"
sid=Tweets_as_A(text)
text=u"レモン!"
sid=Reply_as_B(text, sid)
elif num_now==2 :
text=u"レモンはー！"
sid=Tweets_as_B(text)
text=u"ハチミツ！"
sid=Reply_as_A(text, sid)
text=u"なんでやねん！"
sid=Reply_as_B(text, sid)
#coding:utf-8SCRIPT_TYPE = "InstrumentalityOfMankind" SCRIPT_SHELL = "python" SCRIPT_NAME = "tempFit.py" SCRIPT_HELP_ARG = "" SCRIPT_HELP_SENTENCE = "This is a sample script for" SCRIPT_NUM_ARG = 0 SCRIPT_VERSION = 1.0 ###IMPORT### import sys import os import subprocess ###HELP### if len( sys.argv[1:] ) != SCRIPT_NUM_ARG: print 'Name: ' + SCRIPT_NAME print 'Arguments: ' + SCRIPT_HELP_ARG print 'Explanation: ' + SCRIPT_HELP_SENTENCE sys.exit() ###MAIN### range_min = 2.3 range_max = 10.0 Basic_Model = 'Model_Basic.xcm' Command_error = 'error max 1000 1-29 \n' #デバック用 #Command_error = '\n' for Reg in 'abcd': ###Input file### filename = 'sp.Position_' + Reg + '.FI.rebin.pi' Command_file = \ 'data 1 ' + filename + '\n ' + \ '@' + Basic_Model + '\n' ###Temperature### for Temp in ['1Temp', '2Temp']: if Temp == '1Temp': Command_Temperature = """ newpar 6 1 -1 \n newpar 8 0 -1 \n newpar 9 0 -1 \n newpar 10 0 -1 \n """ Abunds = ['1Solar', 'Free'] Abss = ['Same'] elif Temp == '2Temp': Command_Temperature = """ thaw 6 \n newpar 1 \n thaw 8 \n thaw 9 \n thaw 10 \n """ Abunds = ['1Solar', 'Free', 'Free_Respectively'] Abss = ['Same', 'Free_Respectively'] ###Abund### Command_Abund1 = 'newpar 7=2 \n' for Abund in Abunds: if Abund == '1Solar': Command_Abund = Command_Abund1 + ' newpar 2 1 -1 \n' elif Abund == 'Free': Command_Abund = Command_Abund1 + ' thaw 2 \n' elif Abund == 'Free_Respectively': Command_Abund = Command_Abund1 + """ untie 2,7 \n thaw 2 \n thaw 7 \n """ ###Absorption### for Abs in Abss: if Abs == 'Same': Command_Abs = 'newpar 10 = 5 \n' elif Abs == 'Free_Respectively': Command_Abs = """ thaw 5 \n thaw 10 \n """ ###Output files### output_id = 'Region_' + Reg + '_' + Temp + '_Abund_' + Abund + '_Abs_' + Abs logfile = output_id + '.log' xcmfile = output_id + '.xcm' if os.path.exists( logfile ): os.remove( logfile ) if os.path.exists( xcmfile ): os.remove( xcmfile ) ###Fit### Command_Fit = \ 'ignore **-' + str( range_min ) + '\n ' + \ 'ignore ' + str( range_max ) + '-** \n' + \ """ query yes \n renorm \n fit \n """ + \ Command_error + \ 'log ' + logfile + '\n' + \ 'show \n' + \ Command_error + \ 'log none \n' + \ 'save all ' + xcmfile + ' \n' xsp = subprocess.Popen( 'xspec', stdin = subprocess.PIPE ) xsp.stdin.write( Command_file ) xsp.stdin.write( Command_Temperature ) xsp.stdin.write( Command_Abund ) xsp.stdin.write( Command_Abs ) xsp.stdin.write( Command_Fit ) xsp.stdin.write( 'quit \n \n' ) xsp.wait()
#!/Library/Frameworks/EPD64.framework/Versions/Current/bin/python
SCRIPT_TYPE="InstrumentalityOfMankind"
SCRIPT_SHELL="python"
SCRIPT_NAME="Fit_infile_.py"
SCRIPT_HELP_ARG=""
SCRIPT_HELP_SENTENCE=""
SCRIPT_NUM_ARG=0
SCRIPT_VERSION=1.0
###IMPORT###
import SMM
from scipy.optimize import leastsq
from numpy import *
import sys
###HELP###
if len(sys.argv[1:])!=SCRIPT_NUM_ARG:
print 'Name: '+SCRIPT_NAME
print 'Arguments: '+SCRIPT_HELP_ARG
print 'Explanation: '+SCRIPT_HELP_SENTENCE
sys.exit()
###MAIN###
def SMM_norm_free(norm_nsc,norm_nsd,norm_gb,norm_gd):
return lambda l,b: norm_nsc * SMM.SMM(l,b,'NSC') + norm_nsd * SMM.SMM(l,b,'NSD') + norm_gb * SMM.SMM(l,b,'GB') + norm_gd * SMM.SMM(l,b,'GD')
def SMM_norm_free_array(param,l_array,b_array):
res=[]
num=0
for l in l_array :
b=b_array[num]
res.append(SMM_norm_free(*param)(l,b))
num+=1
return res
def fit_SMM_norm_free(param,l_array,b_array,sb_array,sb_err_array):
errorfunction = lambda p: (array(SMM_norm_free_array(p,l_array,b_array))-array(sb_array))/sb_err_array
p=leastsq(errorfunction,param,full_output=True)
return p
def demo(para=[1.0,2.0,3.0,4.0],para0=[1.1,1.2,1.3,1.4]) :
l_in=linspace(0.0,5.0,10)
b_in=linspace(0.0,5.0,10)
data = SMM_norm_free_array(para,l_in, b_in)
data_err=0.001*array(data)
param=fit_SMM_norm_free(para0,l_in,b_in,data,data_err)
val=param[0]
covar=param[1]
val_err=[]
for num in range(0,4):
val_err.append((covar[num][num])**0.5)
print val,val_err
def fit_SMM_indata(infile):
para0=(1.0,1.0,1.0,1.0)
data=loadtxt(infile)
l_in=data[:,0]
b_in=data[:,1]
sb_in=data[:,2]
sb_in_err=data[:,3]
param=fit_SMM_norm_free(para0,l_in,b_in,sb_in,sb_in_err)
val=param[0]
covar=param[1]
val_err=[]
for num in range(0,4):
val_err.append((covar[num][num])**0.5)
print val,val_err
fit_SMM_indata(str(sys.argv[1]))
#coding:sjis ''' Created on 2011/03/18 @author: uchiyama ''' """ importの使い方・考え方は以下が分かりやすい。 http://1-byte.jp/2010/09/29/learning_python5/ そもそもこのシリーズは図解が多くて良さげなので後で読む。 """ from math import cos,sin,pi, atan, exp, sqrt from scipy import integrate """ 普遍であるべき観測量 ここにグローバルな変数として書くのは正しいのだろうか。 """ #postion of Sgr A* in the Galactic coordinate (degree) L_GC=-0.056 B_GC=-0.046 # 1 pc in the unit of cm #CM_PER_PC=3.09e18 #Position of the Earth from the plane (pc) Z0=16.0 #Distance to Sgr A* (pc) R0=8.5e3 #The rotation angle BD1=15.0 BD2=1.0 """ 座標変換に必要な関数。 本当はscipyの座標変換使えばもっと賢くできる気がする。 """ def sin_d(degree): return sin(2*pi*degree/360.0); def cos_d(degree): return cos(2*pi*degree/360.0); def dl(l): return l-L_GC def db_bulge(b): return b-B_GC-atan(Z0/R0) def db(b): return b-B_GC def r_f(s,l,b): l0=dl(l) b0=db(b) (R0**2+(s*cos_d(b0))**2-2*R0*s*cos_d(l0)*cos_d(b0))**0.5 def z_f(s, l, b): b0=db(b) return s*sin_d(b0) def x_f_bulge(s, l, b): l0=dl(l) b0=db(b) return R0-s*cos_d(b0)*cos_d(l0) def y_f_bulge(s, l, b): l0=dl(l) b0=db_bulge(b) s*cos_d(b0)*sin_d(l0) def z_f_bulge(s, l, b): b0=db_bulge(b) return s*sin_d(b0)+Z0 def X_f(x, y, z): return x*cos_d(BD1)+y*sin_d(BD1) def Y_f(x, y, z): return -1*x*sin_d(BD1)*cos_d(BD2)+y*cos_d(BD1)*cos_d(BD2)+z*sin_d(BD2); def Z_f(x, y, z): return -1*x*sin_d(BD1)*sin_d(BD2)+y*cos_d(BD1)*sin_d(BD2)+z*cos_d(BD2); """ 以下星質量モデル。 基本的にMuno+06に基づく…が数式に誤りたくさん+説明不足。 http://adsabs.harvard.edu/abs/2006ApJS..165..173M Kent+91, Launhardt+02を元に情報を補足、整理。 http://adsabs.harvard.edu/abs/1991ApJ...378..131K http://adsabs.harvard.edu/abs/2002A%26A...384..112L 更に詳しくは http://www-utheal.phys.s.u-tokyo.ac.jp/~uchiyama/wordpress/wp-content/uploads/2011/05/SMM.pdf """ def NSC(x,l,b): norm1=3.3e6 # M_solar pc^{-3} norm2=norm1*8.98e7/3.3e6 n1=2.0 n2=3.0 R_c=0.22 # pc R_lim=6.0 R_cutoff=200.0 r=r_f(x,l,b) z=z_f(x,l,b) """ powや**は（defaultで）使えるのにsqrtはmathをimportしないと使えない。 意味がよく分からないが、そもそもsqrt不要と言われるとそれが正しい気もする。 """ R=sqrt(pow(r,2)+pow(z,2)) """ 元々は3項演算子を使っていたがpythonの3項演算子に慣れずネストの仕方が分からないのでif文で書く。 """ if R < R_lim : return norm1/(1+pow((R/R_c),n1)) elif R < R_cutoff : return norm2/(1+pow((R/R_c),n2)) else: return 0 def NSD(x,l,b): n1=-0.1 n2=-3.5 n3=-10.0 r_1=120.0 r_2=220.0 r_3=2000.0 r_d=1.0 z_d=45.0 norm1=300.0 norm2=norm1*(r_1/r_d)**(-n2+n1) norm3=norm2*(r_2/r_d)**(-n3+n2) r=r_f(x,l,b)+1e-30 z=z_f(x,l,b) if r <=0 : return 0.0 elif r < r_1 : return norm1*(abs(r/r_d)**n1)*exp(-1.0*z/z_d) elif r < r_2 : return norm2*(abs(r/r_d)**n2)*exp(-1.0*z/z_d) elif r < r_3 : return norm3*(abs(r/r_d)**n3)*exp(-1.0*z/z_d) else: return 0.0 def GB(x,l,b): """ Unit of value returned by this function M_solar pc^{-3} """ #norm Unit: M_solar pc^{-3} norm=8.0 #a_x,a_y,a_z : pc a_x=1100.0 a_y=360.0 a_z=220.0 c_para=3.2 c_perp=1.6 x_b=x_f_bulge(x, l, b); y=y_f_bulge(x, l, b); z=z_f_bulge(x, l, b); """ Pythonはべき乘を表す演算子**もあるけどpowもある。 そして小数乘もできるし、負の数を小数乗して複素数を 返すこともできる。 """ R_perp=((abs(X_f(x_b,y,z))/a_x)**c_perp+(abs(Y_f(x_b,y,z))/a_y)**c_perp)**(1/c_perp) Rs=(pow(abs(Z_f(x_b,y,z))/a_z,c_para)+pow(R_perp,c_para))**(1/c_para); return norm*exp(-1*Rs); def GD(x,l,b): """ Unit of this function M_solar pc^{-3} """ #norm Unit: M_solar pc^{-3} norm=5.0 #r_d, z_d : pc r_d=2700.0; z_d=200.0; r=r_f(x, l, b) z=z_f(x, l, b) return norm*exp(-1.0*r/r_d)*exp(-1.0*(z/z_d)) def Stellar_density_f(x, l, b): """ Unit of this function M_solar pc^{-3} """ return NSC(x,l,b)+NSD(x,l,b)+GB(x,l,b)+GD(x,l,b) def SMM(l, b, s=16e3 , region='SMM'): """ Pythonで関数はオブジェクトなので下記の様なことができてしまう。 正直キモイ。ポインタを渡したくなる…。以下を参照。 http://python.g.hatena.ne.jp/muscovyduck/20080707/p2 """ if region=='NSC': func=NSC elif region=='NSD': func=NSD elif region=='GD': func=GD elif region=='GB': func=GB elif region=='SMM': func=SMM else : exit("3rd argument must be one of NSC/NSD/GD/GB/SMM.") """ SciPy(さいぱい。すしぱいと読んでいた…)を使った積分。 http://docs.scipy.org/doc/scipy/scipy-ref.pdf の1.4.1を見る。 lambda式の意図は http://atkonn.blogspot.com/2008/02/python-python27-lambda.html をみるとなんとなく分かる。 """ res=integrate.quad(lambda x: func(x,l,b), 0.0, s)[0] """ The unit of returned value is M_solar/pc^2/str """ return res/(4.0*pi);
<pre class="prettyprint" id="python"> #coding:utf-8 import pylab print "Hello, world!" x=pylab.arange(0.0,1.0,0.01) y=pylab.sin(2*pylab.pi*x) figplot=pylab.subplot(111) figplot.plot(x,y) pylab.show() </pre>んで
#coding:utf-8 import pylab print "Hello, world!" x=pylab.arange(0.0,1.0,0.01) y=pylab.sin(2*pylab.pi*x) figplot=pylab.subplot(111) figplot.plot(x,y) pylab.show()
#coding:utf-8
"""
Gnuplot.pyが何故かeclipseで補完されない。
ライブラリに手動で設定したのだが。
何故か gnuplot_Suites が出てくる。
"""
import Gnuplot
import numpy as np
in_file='Data_all.txt'
in_txt_name='Data_GalacticPlane.txt'
out_file='Gnuplot_test.ps'
out_file2='Gnuplot_test2.ps'
"""
Pythonらしい?方法でファイル読み込み
"""
input_file=open(in_file,"r")
data_list=[]
Fe=[]
l=[]
line_data_cols=[4,7,10]
line_ids=[0,3,6]
for id in range(9):
Fe.append([])
for input_line in input_file:
if input_line.strip().split(',')[14]==' l':
l.append(float(input_line.strip().split(',')[2]))
for line_col_num in line_data_cols:
"""
line_id:
0-2 Fe I
3-5 Fe XXV
6-8 Fe XXVI
"""
line_id=line_col_num-4
value=float(input_line.strip().split(',')[line_col_num])
lower=float(input_line.strip().split(',')[line_col_num+1])
upper=float(input_line.strip().split(',')[line_col_num+2])
Fe[line_id].append(value*1e-7)
Fe[line_id+1].append(lower*1e-7)
Fe[line_id+2].append(upper*1e-7)
"""
numpyのloadtxtを使ってもう少し賢くやる
http://d.hatena.ne.jp/y_n_c/20091117/1258423212
を参考にする。
便利だが今回のようにデータを各行毎にチェックした上でplotしなくては
ならない場合は向かないのかもしれない。
l=-0.05のデータのみ抜き出したファイルを準備しておく。
"""
l_txt=np.loadtxt(in_txt_name, delimiter=',', usecols=[3],unpack=True)
"""
(np.)arrayとlistのどちらを使うべきかでよく混乱するが「自分の印象としては」
arrayは静的な数字そのものでlistは動的な変数。
・arrayにはappendはないけど、sin(array)みたいなことをやって、更にarrayを作れる。
・array*2=array[2*a,2*b]だが[a,b]*2=[a,b,a,b]
・そもそもarrayには数字以外は入らない。
"""
Fe_txt=[]
Fe_max_txt=[]
Fe_min_txt=[]
for id in range(3):
Fe_txt.append([])
Fe_max_txt.append([])
Fe_min_txt.append([])
l_txt=np.loadtxt(in_txt_name, delimiter=',', usecols=[2],unpack=True)
for num in range(3):
Fe_txt[num],Fe_min_txt[num],Fe_max_txt[num]=np.loadtxt(in_txt_name, delimiter=',', usecols=[3*num+4,3*num+5,3*num+6],unpack=True)
"""
Gnuplot.pyでプロット。
・Gnuplot.pyの正式なマニュアル
http://gnuplot-py.sourceforge.net/doc/
http://gnuplot-py.sourceforge.net/doc/Gnuplot/_Gnuplot/Gnuplot.html
http://www-acc.kek.jp/EPICS_Gr/products/gnuplot/Gnuplot_py.pdf
が日本語で詳しくて有り難い。
・withはpythonでは予約語なのでwith_を使う。
・Gnuplot.pyを使ってfittingも出来なくは無さそうだが、
今一つ頭が悪い感じ。少なくともpython内で定義した関数で
fitは出来そうにない。
http://osdir.com/ml/python.ipython.user/2003-07/msg00009.html
"""
"""
Python I/Oで読んだデータをプロット。
"""
gp1=Gnuplot.Gnuplot()
lw='3'
ps='1.5'
data1=[]
for num in range(3):
data1.append([])
for num in range(3):
if num==0:
line_name='Fe I K{/Symbol a}'
pt='7'
lt='0'
elif num==1:
line_name='Fe XXV K{/Symbol a}'
pt='6'
lt='1'
elif num==2:
line_name='Fe XXVI K{/Symbol a}'
pt='9'
lt='3'
"""
何故か titleだけ with_でなくてtitleでGnuplot.Dataに渡す。
"""
with_option='e lw '+lw+' pt '+pt
data1[num]=Gnuplot.Data(l,Fe[3*num],Fe[3*num+1],Fe[3*num+2],with_=with_option,title=line_name)
gp1('set key')
gp1('set log y')
gp1('set te unknown')
gp1.xlabel('Galactic longitude (degree)')
gp1.ylabel('Intensity (ph^{-1} cm^{-2} arcmin^{-2} s^{-1})')
gp1.title('Fe K{/Symbol a} profile along the Galactic plane')
gp1.plot(data1[0], xrange='[2.2:-3.2]', yrange='[1e-8:1e-5]')
for num in [1,2]:
gp1.replot(data1[num])
gp1('set format y "10^{%L}"')
gp1('set te po co solid enhanced "" 24')
gp1('set ou '+'"'+out_file+'"')
gp1.replot()
gp1('set ou')
"""
loadtxtで読んだデータをプロット。
"""
data2=[]
for num in range(3):
data2.append([])
for num in range(3):
if num==0:
line_name='Fe I K{/Symbol a}'
pt='7'
lt='0'
elif num==1:
line_name='Fe XXV K{/Symbol a}'
pt='6'
lt='1'
elif num==2:
line_name='Fe XXVI K{/Symbol a}'
pt='9'
lt='3'
with_option='e lw '+lw+' pt '+pt
data2[num]=Gnuplot.Data(l_txt,Fe_txt[num],Fe_min_txt[num],Fe_max_txt[num],with_=with_option,title=line_name)
gp1('set te unknown')
gp1.plot(data2[0], xrange='[2.2:-3.2]', yrange='[1e-1:1e2]')
for num in [1,2]:
gp1.replot(data2[num])
gp1.ylabel('Intensity (10^{-7} ph^{-1} cm^{-2} arcmin^{-2} s^{-1})')
gp1('set te po mono dashed enhanced "" 24')
gp1('set ou '+'"'+out_file2+'"')
gp1.replot()
gp1('set ou')
#coding:utf-8
import pylab
input_file=open("../../Data_all.txt","r")
data_list=[]
Fe=[]
l=[]
line_data_cols=[4,7,10]
line_ids=[0,3,6]
for id in range(9):
Fe.append([])
for input_line in input_file:
if input_line.strip().split(',')[14]==' l':
l.append(float(input_line.strip().split(',')[2]))
for line_col_num in line_data_cols:
"""
line_id:
0-2 Fe I
3-5 Fe XXV
6-8 Fe XXVI
"""
line_id=line_col_num-4
value=float(input_line.strip().split(',')[line_col_num])
lower=float(input_line.strip().split(',')[line_col_num+1])
upper=float(input_line.strip().split(',')[line_col_num+2])
Fe[line_id].append(value*1e-7)
Fe[line_id+1].append((value-lower)*1e-7)
Fe[line_id+2].append((upper-value)*1e-7)
"""
Fontの設定の仕方は
http://www.scipy.org/Cookbook/Matplotlib/LaTeX_Examples
(rcParams.updateの使い方)
具体的にパラメータとして何が使えるかは
http://matplotlib.sourceforge.net/users/customizing.html#a-sample-matplotlibrc-file
とかを見る。
"""
fparams = {'backend': 'ps',
'axes.labelsize': 20,
'suptitle.fontsize': 20,
'legend.fontsize': 20,
'xtick.labelsize': 18,
'ytick.labelsize': 18,
'text.usetex': True,
'font.family': 'serif',
'font.serif': 'Times New Roman',
'text.usetex' : True}
pylab.rcParams.update(fparams)
l_dist=pylab.subplot(111)
"""
ログスケールの図の作り方は下記。
http://matplotlib.sourceforge.net/examples/pylab_examples/log_demo.html#pylab-examples-example-code-log-demo-py
"""
l_dist.set_yscale("log", nonposy='clip')
l_dist.set_ylim(ymin=1e-8,ymax=1e-5)
l_dist.set_xlim(xmin=2.2,xmax=-3.3)
"""
TeXの書式を使う方法は下記。
http://www.scipy.org/Cookbook/Matplotlib/UsingTex
"""
for line_id in line_ids:
if line_id==0:
title=r'Fe I K$\alpha$'
elif line_id==3:
title=r'Fe XXV K$\alpha$'
elif line_id==6:
title=r'Fe XXVI K$\alpha$'
"""
エラーのついた図の作り方は下記。
http://matplotlib.sourceforge.net/api/pyplot_api.html#matplotlib.pyplot.errorbar
"""
l_dist.errorbar(l,Fe[line_id], yerr=[Fe[line_id+1],Fe[line_id+2]], fmt='.', label=title, elinewidth=3 )
"""
凡例(Gnuplotで言うところのKey)の作り方は下記。
http://matplotlib.sourceforge.net/api/artist_api.html#matplotlib.legend.Legend
handletextpadは凡例中のマークと文章の間のスペース
"""
l_dist.legend(numpoints=1,frameon=False,markerscale=1.5,handletextpad=-0.5)
"""
suptitleのみ fontsizeが必要な理由が自分には意味不明。
"""
pylab.suptitle(r'Profile of the Fe K$\alpha$ lines along the Galactic longitude ($b=-0^{\circ}.05$)',size=20)
pylab.xlabel(r'Galactic longitude $l$ (degree)')
pylab.ylabel(r'Intensity (ph s$^{-1}$ cm$^{-2}$ arcmin$^{-2}$ )')
"""
表示した図を消さないためのpause
"""
raw_input()
pylab.savefig('GCXE.eps', format='eps')
