import ddf.minim.*; import ddf.minim.ugens.*; Minim minim; AudioOutput out; Oscil sin1, sin2; float f=440; float A=1.0; float trig=0.5; void setup() { size(512, 200); smooth(); minim = new Minim(this); out = minim.getLineOut(Minim.MONO, 1024); sin1=new Oscil(f, A, Waves.SINE); sin2=new Oscil(f*2, A*0.5, Waves.SINE); sin1.setPhase(0); sin2.setPhase(0); sin1.patch(out); sin2.patch(out); } void draw() { background(255); stroke(0); float x=0, Amp_old=999; int flag=0; for (int i = 0; i < out.bufferSize()-1; i++) { float Amp=out.left.get(i); if ((Amp_old=trig)||flag==1) { point(x, height/2 + Amp*height/4); flag=1; x+=2; } Amp_old=Amp; } } void mouseMoved() { float phase=map(mouseX, 0, width, 0, 1); sin2.setPhase(phase); } void stop() { out.close(); minim.stop(); super.stop(); }
from pyorbital import orbital ＃これを追加すると動く（外すと動かない）のは何故
from pyorbital import orbital
import pyorbital
import numpy
import datetime
#Sat info
satname='ITUPSAT 1'
f=437.325
tlef='TLE.txt'
#Position Info
#静岡大学
lon=138.431704
lat=34.965720
alt=90
#Obs Info
start_jst=datetime.datetime(2018,02,16,15,25,00)
end_jst=datetime.datetime(2018,02,16,15,45,00)
#Conf
c=299792.458 # km/s
dt=datetime.timedelta(0,30)
def jst2utc(time) :
return time-datetime.timedelta(0,9.*60*60,0)
def obsf_ICR6(obsf) :
f0=obsf*1000.
f1=int(obsf*100.)*10
df=f0-f1
if df <2.5 :
return f1/1000.
elif df <7.5 :
return (f1+5.)/1000.
else :
return (f1+10.)/1000.
tle=pyorbital.tlefile.read(satname,tlef)
orb=pyorbital.orbital.Orbital(satname,tlef)
print "TLE Epoch:",",",str(datetime.datetime(2000+int(tle.epoch_year),1,1)+datetime.timedelta(days=tle.epoch_day)).split(' ')[0]
print "Date:",",",str(start_jst).split(' ')[0]
print "Lon:",",",lon
print "Lat:",",",lat
print "Alt:",",",alt
print ""
jst=start_jst
while jst < end_jst :
now=jst2utc(jst)
obs_pos=numpy.array(pyorbital.astronomy.observer_position(now,lon,lat,alt)[0])
sat_pos,sat_vel=orb.get_position(now,False)
range_dis=numpy.linalg.norm(sat_pos-obs_pos)
range_vec=(sat_pos-obs_pos)/range_dis
rang_rate=numpy.dot(sat_vel,range_vec)
Az,El=orb.get_observer_look(now,lon,lat,alt)
obsf=f*(1-rang_rate/c)
print str(jst).split(' ')[1],",",Az,",",El,",",obsf_ICR6(obsf),",",obsf
jst+=dt
int a,b,c; float n1=0; float n2=0; void setup(){ size(1000,20); } void draw(){ a=int(random(1,7)); b=int(random(1,7)); c=int(random(1,7)); if(log(a+b)/log(0.25)>log(c)/log(0.5)) n1+=1; if((pow(2,a)+pow(2,b)+pow(2,c))%3==0) n2+=1; background(255); fill(0); text(str(frameCount),width*0.1,height/2); text(str(n1/frameCount),width*0.5,height/2); text(str(n2/frameCount),width*0.8,height/2); }数えた。コンピュータ便利。
int a, b, c; float n1=0; float n2=0; float n=0; size(1000, 20); for (a=1; a<7; a+=1) { for (b=1; b<7; b+=1) { for (c=1; c<7; c+=1) { if (log(a+b)/log(0.25)>log(c)/log(0.5)) n1+=1; if ((pow(2, a)+pow(2, b)+pow(2, c))%3==0) n2+=1; n+=1; } } } background(255); fill(0); text(str(n), width*0.1, height/2); text(str(n1/n), width*0.5, height/2); text(str(n2/n), width*0.8, height/2);
PImage myPhoto; int myPhotoWidth, myPhotoHeight; void setup() { size(500,750); //Thanks for http://gahag.net/004881-high-school-girl/ myPhoto = loadImage("Test.jpg"); myPhotoWidth = myPhoto.width; myPhotoHeight = myPhoto.height; noStroke(); smooth(); background(0); myPhoto.loadPixels(); } void draw() { background(0,0,100); int x; int y; for(x=0;x<myPhoto.width;x+=5){ for(y=0;y<myPhoto.height;y+=5){ strokeWeight(random(1,3)); float f=random(0.5,1); stroke(red(myPhoto.pixels[y*width + x])*f,green(myPhoto.pixels[y*width + x])*f,blue(myPhoto.pixels[y*width + x])*f); float r0=random(5, 15); float r1=random(5, 15); float th0=random(0, 2*PI); float th1=random(0, 2*PI); float dx=random(-5, 5); float dy=random(-5, 5); line(r0*cos(th0)+x+dx, r0*sin(th0)+y+dy, r1*cos(th1)+x+dx, r1*sin(th1)+y+dy); }} }
PImage img;
int R=200, G=200, B=200; void setup() { size( 800, 800 ); background( R, G, B ); //Thanks for Yume Yume Iro TOWN http://www.poipoi.com/yakko/cgi-bin/sb/log/eid1737.html img = loadImage("siruetto20080205a.png"); stroke(255, 0.0); noCursor(); } void draw() { background(R, G, B); float x=mouseX-280; float y=mouseY-500; for (int i=0; i<10; i++) { fill(R, G, B, 30); noStroke(); rect(0, 0, width, height); float xi=map(i, 0, 10, width/2, x); float yi=map(i, 0, 10, height/2, y); float f=map(i, 0, 10, 0.1, 1); image(img, xi, yi, 560*f, 1000*f); } } void mouseMoved() { int a=int(random(1, 4)); int pm; if (R+G+B<50) { pm=1; } else if (random(0, 1.0)<0.5) { pm=-1; } else { pm=1; } if (a==1) { if (R==200) { R=199; } else if (R==1) { R=2; } else { R=R+pm; } } if (a==2) { if (G==200) { G=199; } else if (G==1) { G=2; } else { G=G+pm; } } if (a==3) { if (B==200) { B=199; } else if (B==1) { B=2; } else { B=B+pm; } } }
PImage img; void setup() { size( 600, 600 ); background( 0, 0, 255 ); //Thanks for Yume Yume Iro TOWN http://www.poipoi.com/yakko/cgi-bin/sb/log/eid1737.html img = loadImage("siruetto20080205a.png"); stroke(255, 0.0); } void draw() { background(0, 0, 255); image(img, mouseX-140, mouseY-250, 280, 500); //image(img, mouseX-120, mouseY-80, 280, 500); beginShape(); for (int i=0; i<300; i++) { stroke(random(0, 255), 0, 0); strokeWeight( random(1, 3) ); float r0=random(0, 50); float r1=random(0, 50); float th0=random(0, 2*PI); float th1=random(0, 2*PI); line(r0*cos(th0)+mouseX-25, r0*sin(th0)+mouseY-170, r1*cos(th1)+mouseX-25, r1*sin(th1)+mouseY-170); } }
PImage img; void setup() { size( 800, 800 ); background( 255 ); //Thanks for Yume Yume Iro TOWN http://www.poipoi.com/yakko/cgi-bin/sb/log/eid1737.html img = loadImage("siruetto20080205a.png"); strokeWeight( 6 ); } void draw() { if (frameCount%3==0) { float f=random(0.1, 0.8); image(img, random(0, 500), random(0, 500), 560*f, 1000*f); } stroke( random(255), random(255), random(255) ); fill( random(255), random(255), random(255)); beginShape(); int x0=radomxy(); int y0=radomxy(); vertex( x0, y0 ); bezierVertex( radomxy(), radomxy(), radomxy(), radomxy(), radomxy(), radomxy() ); int i=0; while (i==0) { if (random(6)>3) { i=1; } bezierVertex( radomxy(), radomxy(), radomxy(), radomxy(), radomxy(), radomxy() ); } bezierVertex( radomxy(), radomxy(), radomxy(), radomxy(), x0, y0 ); endShape(); } void mousePressed() { //background( 255 ); delay(1000); } int radomxy() { return int(random(-100, height+100)); }
int x, y, vx, vy; void setup() { smooth(); size(800, 800); x=width/2; y=height/2; background(255); int v=20; vx=int(random(-1*v, v)); int a; if (random(1)>0.5) { a=1; } else { a=-1; } vy=a*int(sqrt(pow(v, 2)-pow(vx, 2))); } void draw() { float a=random(1, 4); int xnew=x+int(a*vx); int ynew=y+int(a*vy); if (xnew<0) { xnew=0; vx=-1*vx; vy=int(vy*random(0.9, 1.1)); } if (xnew>width) { xnew=width; vx=-1*vx; vy=int(vy*random(0.9, 1.1)); } if (ynew<0) { ynew=0; vy=-1*vy; vx=int(vx*random(0.9, 1.1)); } if (ynew>height) { ynew=height; vy=-1*vy; vx=int(vx*random(0.8, 1.2)); } if (vx*vx<30) { vx=int(random(-20, 20)); } if (vy*vy<30) { vy=int(random(-20, 20)); } stroke(random(255), random(255), random(255)); strokeWeight(random(30)); line(x, y, xnew, ynew); x=xnew; y=ynew; if (frameCount%10==0) { fill(255, 20); noStroke(); rect(0, 0, width, height); } }
float theta=0; int n=0; void setup() { size(1000, 1000); background(255); noLoop(); } void draw() { float r=width/3; float re=width/8; float rm=width/10; background(255); stroke(100); fill(0, 0, 255); ellipse(width/2, height/2, re, re); stroke(0); fill(255, 255, 0); ellipse(width/2+r*cos(theta), height/2+r*sin(theta), rm, rm); fill(255, 0, 0); ellipse(width/2+r*cos(theta)+0.35*rm*cos(PI+1*theta), height/2+r*sin(theta)+0.35*rm*sin(PI+1*theta), rm*0.2, rm*0.2); theta+=PI/200; } void mousePressed() { if (n==0) { loop(); n=1; } else { noLoop(); n=0; } }
float n=1; float s=0; void setup() { size(300, 350); background(255); noFill(); stroke(0, 2555, 0); ellipse(0, 0, 600, 600); } void draw() { float x=random(0.0, 1.0); float y=random(0.0, 1.0); if (pow(x, 2)+pow(y, 2)<1.0) { stroke(255, 0, 0); s+=1; } else { stroke(0, 0, 255); } ellipse(x*300, y*300, 1, 1); stroke(255); fill(255); rect(0, 310, 300, 50); fill(0); text(str(int(n)), 5, 330); text(str(s/n*4.0), 150, 330); n++; } void mousePressed() { n=1; s=0; background(255); noFill(); stroke(0, 2555, 0); ellipse(0, 0, 600, 600); }
float n=1; void setup(){ size(300,350); background(255); } void draw(){ float s=0; background(255); fill(255); noFill(); stroke(255,0,0); ellipse(0,0,600,600); float x0=0; float y0=1.0; for(int m=1;m<n+1;m++){ float x=m/n; float y=sqrt(1-pow((x),2.0)); stroke(0); line(x*300,y*300,x0*300,y0*300); s=s+sqrt(pow((x-x0),2.0)+pow((y-y0),2.0)); x0=x; y0=y; } n+=int(n/10+1); if(n<300){ delay(int(1000/sqrt(n))); } fill(0); textSize(15); text(str(int(n)),5,330); text(str(s*2.0),150,330); if(n>10000000){ noLoop(); } } void mousePressed(){ n=1; loop(); }
float n=1; void setup(){ size(300,350); background(255); } void draw(){ float s=0; background(255); fill(255); for(int m=0;m<n;m++){ float x=m/n; float f=sqrt(1-pow((x),2.0)); stroke(0); rect(x*300,0,1/n*300,f*300); noFill(); stroke(255,0,0); ellipse(0,0,600,600); s=s+f*1/n; } n+=int(n/10+1); if(n<300){ delay(int(1000/sqrt(n))); } fill(0); textSize(15); text(str(int(n)),5,330); text(str(s*4.0),150,330); if(n>10000000){ noLoop(); } } void mousePressed(){ n=1; loop(); }
float glevel=0.8;
float x;
float xb, yb, vy=0, vxw;
float headsize;
int t=0;
int ball_start=150;
float g;
float[] xbo=new float[500];
float[] ybo=new float[500];
void setup() {
x=100;
size(1100, 700);
colorMode(RGB, 256);
background(255, 255, 255);
headsize=width*0.03;
frameRate(50);
xb=0;
yb=height*glevel*0.15;
g=0.05;
vxw=width*0.002;
for (int i=0; i<500 i="" p=""> xbo[i]=0;
ybo[i]=0;
}
}
void draw() {
background(255, 255, 255);
//smooth() ;
//World
stroke(0, 0, 0);
strokeWeight(1);
line(0, height*glevel, width, height*glevel);
noStroke();
//point
fill(255, 0, 0);
ellipse(x, height*glevel, width*0.01, width*0.01);
//man
fill(0, 0, 0);
//head
ellipse(x+width*0.1, height*glevel*0.71, headsize, headsize);
//body
rect(x+width*(0.1-0.025), height*glevel*0.75, width*(0.05), height*0.1);
//leg
rect(x+width*(0.1-0.025), height*glevel*0.89, width*(0.02), height*0.085);
rect(x+width*(0.1+0.005), height*glevel*0.89, width*(0.02), height*0.085);
//arm
rect(x+width*(0.1-0.045), height*glevel*0.75, width*(0.015), height*0.08);
rect(x+width*(0.1+0.03), height*glevel*0.75, width*(0.015), height*0.08);
//Bar
strokeWeight(10);
stroke(183, 65, 14);
line(x+width*0.2, height*glevel, x+width*0.2, height*glevel*0.1);
line(x, height*glevel*0.1, x+width*0.2, height*glevel*0.1);
//Ball Shadow
for (int i=1; i if ( i%20==0) {
stroke(150);
strokeWeight(1);
noFill();
ellipse(xbo[i], ybo[i], width*0.03, width*0.03);
}
}
//Ball
xbo[t]=xb;
ybo[t]=yb;
if (t>ball_start & yb<=height*glevel) {
vy=vy+g;
yb=yb+vy;
//下の1文を入れるか入れないかが両者の差。
xb=x;
} else {
xb=x;
}
if (yb>=height*glevel) {
vy=0;
yb=height*glevel;
xb=xb+vxw;
}
noStroke();
fill(0, 0, 255);
ellipse(xb, yb, width*0.03, width*0.03);
x=x+vxw;
print(t, "\n");
if (t<499 p=""> t++;
}
}
import ddf.minim.*;
import ddf.minim.signals.*;
Minim minim;
AudioOutput out;
SineWave[] sine;
float waveH = 50;
int n=30;
void setup()
{
size(512, 200);
smooth();
minim = new Minim(this);
sine=new SineWave[n];
out = minim.getLineOut(Minim.STEREO);
int i;
for(i=0;i < n;i++){
int m=i+1;
sine[i] = new SineWave(440*m, 2*(1-cos(m*PI))/(m*PI*440)*100, out.sampleRate());
}
for(i=0;i < n;i++){
out.addSignal(sine[i]);
}
}
void draw()
{
background(255);
stroke(0);
for(int i = 0; i < out.bufferSize()-1; i++)
{
point(i, 50 + out.left.get(i)*waveH);
point(i, 150 + out.right.get(i)*waveH);
}
}
void mouseMoved() {
int t = int(map(mouseX, 0, width, 1, n));
for(int i=0;i < t;i++){
int m=i+1;
sine[i].setAmp(2*(1-cos(m*PI))/(m*PI*440.)*100);
}
for(int i=t;i < n;i++){
sine[i].setAmp(0);
}
}
void stop()
{
out.close();
minim.stop();
super.stop();
}
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')
