2017年6月26日月曜日

大学入試問題をモンテカルロで解く

情報処理の課題として、「大学入試問題をモンテカルロで解く」というのは面白いかと思い、こちらの問題について作成。しかし、これなら素直に数え上げた方が遥かに速い(汗)。


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);

2017年5月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);
   }}
}

2017年4月30日日曜日

色々統合した練習

もう少し賢く書けそう。こちら
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;
    }
  } 
}

2017年4月29日土曜日

マウスの使い方練習

マウスの使い方練習。IRHH 2
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));
}

2017年4月26日水曜日

なにかとむなしい

こちら

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);
  }
}

2017年4月19日水曜日

月の自転・公転周期デモ

月の自転・公転周期の説明デモ。動いているのはこちら。意図せずもドラoもんカラー。
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;
  }
}  

2017年4月14日金曜日

マル9円周率その3

意味無しモンテカルロ。こちら
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);
}  

マル9な円周率その2

その2。動いているのはこちら
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();
}  

2017年4月13日木曜日

マル9な円周率計算デモ

頭の悪い円周率計算デモ。Processing.jsで動かすとこちら
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();
}  

2017年4月7日金曜日

地球の自転に対するアリストテレス or ニュートン的描像

天文学の授業で「地動説に対する反論」を説明する為にProcessingで作成したアニメーション。こういう教材を簡単に作れるので、本当はどこかの授業で学生さんにProcessingを触ってもらいたいところ。Processing.jsで動かしたアリストテレス or ニュートン的な描像はそれぞれ、こちらこちら
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++;

  }

}

2017年4月5日水曜日

Processingで矩形波のフーリエ変換の可視・音声化デモ

こちらを参考にさせていただきました。学生さんの卒論に向けたProcessingでのフーリエ合成(というかMinimの使い方)練習。Processing、本当に楽しい...。

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();
}

2016年11月26日土曜日

Processingでブラウン運動

実験の授業の参考のためにProcessingでブラウン運動 (というか酔歩) を再現するプログラムを書いたところ、やたらとアーティスティックな感じになって愉快。
こちらを参考にしました。
Processing.jsでweb上で動かすとこんな感じこんな感じ。

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;
}

2016年10月5日水曜日

2HD with Arduno

最悪なサンプル、その2。ArdunoとProcessingの連携。やってみて、学生さんに全部やってもらうのはしんどそう、と分かったのは収穫ではあるが...。
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);
 }
}

2016年10月3日月曜日

2HD

最悪なサンプル、その1。 学生さんとProcessingをやっていると、つい楽しくカッとなって作ってしまった...。 サンプルのPerspectiveほぼそのまま。
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);
  }
}

2016年5月23日月曜日

Bashスクリプトで中心極限定理

中心極限定理実験 Bashスクリプト。学生さんの(スクリプトの書き方も含めた) 勉強用。
# 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
###

2016年5月10日火曜日

Labviewでジャンケン

ついカッとなって作った。

https://www.evernote.com/l/ACwx-FddrV5AoIIiAnkGQwOoR1iwTALO-3k
まあ、学生さんへのLaboviewのサンプルということで (汗).




http://rightway7049.com/wp-content/uploads/2015/10/52128c1d.gif
http://dic.nicovideo.jp/oekaki/509304.png
http://www.geocities.jp/kankei33/2006010210_1113818290.jpg
http://getnews.jp/img/archives/imp/and_136717.jpg
http://file.onepiece.ria10.com/IMG_0463.JPG
http://s.eximg.jp/exnews/feed/Kotaku/Kotaku_201005_sf2_doraemon.jpg

2016年4月19日火曜日

お察し下さい

ちくしょう...
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()))

2016年4月7日木曜日

Bashで複数行をコメントアウトする

#! /bin/bash

echo "Inu"
echo "Saru"
echo "Kiji"

:<<"_Test_"
echo "Neko"
ls

xspec<<EOF
exit
EOF
_TEST_

2016年3月31日木曜日

Arf_BG_corrector.py

お察し下さい (pyfits練習) 第2弾。前よりは多少ファイルI/Oの書き方が賢くなっている...かも。あくまで個人的に作ったスクリプト。
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()

2016年3月27日日曜日

openpyxlを使ってPythonでExcelファイルを読み書きする練習

openpyxlを使ってPythonでExcelファイルを読み書きする練習。
###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()

2016年2月24日水曜日

PI Shifter for Suzaku XIS

お察し下さい。個人的に作ったスクリプト。pyfitsの練習というか。

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="Modiefy 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()

2015年10月18日日曜日

Pythonで文字列と数値を含んだテキストファイルの読み込み

こんなファイルを読みたい。 本当はnumpyのloadtxtでdtypeを指定する (or 読みたい列だけ読む) のが賢い気は果てしなくする。
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 適材適所。

2015年7月7日火曜日

N人のクラスで同じ誕生日の人がM人以上いる確率

学生さんのレポートでの質問に答えるために作った。40人クラスで3人以上同じ誕生日の確率は6.7%くらい。
###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()

2015年7月1日水曜日

篠本先生のじゃんけんマシン in Python

篠本先生の教科書で昔勉強したパーセプトロンのじゃんけんマシン (ソースはここを参照しました) をPythonで書き直してみた。 もし情報処理の授業を持つ事になったら、Pythonで機械学習とか面白いかと思ったのだが、ハードルが高い (&篠本先生の教科書は既に入手困難) か...。 GUIをつけてみたいところ。
###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()

2015年1月8日木曜日

ついカッとして作った

今は反省しています...。...辞書の使い方の練習、とか(汗)。
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()

2014年12月24日水曜日

ヘルスライスでバーサーク

『カザンの闘技場』でもらえるヘルスライス(42D)でバーサークするとどの位ヒットが出るのか?」という子供の頃の疑問に答えるスクリプトを誘惑に負けて作成。関数を再帰的に使うことの練習ということで(汗)。

100万回トライして
平均 :  357
標準偏差 :  23
最大 :  501
最小 :  268

非バーサーク時(期待値 147)の2.4倍。
分布でなんか変なヒゲが生えているのは何なのだろう...?



ハイパーバーサークだと1万回トライして
平均 :  499
標準偏差 :  51
最大 :  765
最小 :  348

...仕事します。

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()

2014年12月23日火曜日

中心極限定理をテストするスクリプト

学生さんに宿題で出したので自分でも書いてみた。
Python (matplotlib, numpy) 便利。
Enthought Canopyも便利。

# -*- 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()

 

2013年7月1日月曜日

Twitterで勝手に掛け合いをするBot的な何か

ついカッとなって、週末作った。
今は反省している。 (これをcrontabで定期的に走らせたりすると・・・ゴホゴホ)
勢いでやったからよく覚えていないがこれを使ったっぽい。

#!/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)

2012年12月26日水曜日

LabViewで拡大

は、できないそうだがMacだとcontrol+マウスホイールで幸せになれる。

2012年4月28日土曜日

PythonでXspec

Shellスクリプトからの卒業?を目指してPythonでXSPECの解析用スクリプトを書いてみる。
サワDに教えてもらった方法を使用。

やりたいことは「たくさんのRegionのスペクトルをビミョーに条件を変えながらモデルFit、エラーを計算」。読み込んでいるModel_Basic.xcmはこんな感じ。

やってみた感じとしては、(便利なところも多いが)今回の様に「たくさんのスペクトルをモデルFit、エラー計算」するだけなら、Bashスクリプトで十分(かつ便利)な気がする。
numpyとかscipyとか連携させたい時は良いけど、xspecを操作する部分までpythonでやる必要があるかは謎。

xsp.stdin.write( Command_Fit )とかは、もはや「お約束」なので、関数&モジュール化して再利用しやすくしておくと楽になるのかもしれない。

PyXspecを使えば、別次元の扉が開ける気もするが、import するとFatal Python error・・・何故だ・・・。

他の人がどうやっているのか & sherpaを使った場合を知りたいところ。


#coding:utf-8
SCRIPT_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()

2012年4月27日金曜日

WebExが"Connecting..."で停止する

TV会議で使っているCiscoのWebExが"Connecting..."で停止してしまって困っていたが ~/Library/Application Support/WebEx Folderを丸ごと削除すると問題解決。

2012年4月25日水曜日

set table

新しいgnuplotで"set terminal table"が使えなくて困っていたが"set table"で同じことができるようだ。

http://www.gnuplot.info/docs_4.2/gnuplot.html#x1-27000043.61

2011年8月31日水曜日

立体角

立体角をなかなか理解してもらえなかったので見つけてきたサイトとか。


----

立体角の比較的分かりやすい定義の説明です。

http://www.earliness-si-unit.com/hozyo_02.html

以下、電磁気の講義ノートetc.に現れる立体角の説明です。
いずれも分り易いと思うのでよく読んで下さい。
立体角はガウスの定理を習った時(マクスウェル方程式からクーロンの法則を導出した時)
に「必ず」習っている筈の基本概念なのでよく理解しておいて下さい。

http://www.phys.u-ryukyu.ac.jp/~maeno/cgi-bin/pukiwiki/index.php?%C5%C5%BC%A7%B5%A4%B3%D82007%C7%AF%C5%D9%C2%E84%B2%F3#q9a03723

http://www.erp.sophia.ac.jp/Projects/ocw/faculty/st/stCl/common_st/10SSSCT107E1/10012SCT107E1.pdf

http://hooktail.sub.jp/vectoranalysis/GaussSolidAngle/

2011年7月27日水曜日

2次元 SMM fit

#!/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]))

2011年7月25日月曜日

各XIS (0,1,3) とHXDのNormが食い違う理由

XISとHXDのNormが食い違う理由についての考察?のようなもの。
後輩から質問を受けたのでまとめた。
間違えているかもしれないので注意。指摘していただけると助かります。


-------

以下、簡単に各XIS (0,1,3) とHXDのNormが食い違う理由になりうる較正上の
不定性について簡単に書いておきます。

1. XISのnon-X-ray backgroundの不定性
宇宙線由来のバックグラウンドは観測(衛星軌道)毎・検出器毎、更に検出器上の
場所毎に異なります。XIS 0,1,3間のnormalizationが観測毎に異なるのは主にこのNXBの
揺らぎによる筈です。暗い天体のハードな帯域 (5-10 keVあたり) のfluxはこの不
定性の影響を大きく受けます。 BGを視野内の暗い所から引くという解析では、
この不定性の影響が大きくなります。 xisnxbgenというftoolを使い、「(1) SRCとBG
の領域毎にNXBのスペクトルを作り、(2) SRCとBGの各々NXBを差し引いた上で更に
SRCからBGを引く」といったことを行えば、NXBの不定性は小さくなりますが少し
統計誤差が大きくなります。

2. XISのcontaminationの較正の不定性
XISには可視光遮断フィルターがついていますが、この表面にX線を吸収してしまう汚染
物質(ゴム等からアウトガス)が付着しつつあることが知られています。この汚染物質は
観測時期が進むにつれ増加しており、検出器毎、更に検出器上の場所毎に異なります。
この汚染物質量の較正の不定性はXIS 0,1,3間のnormalizationが異なる原因になります。
ソフトな帯域 (0.5-2 keVあたり) のfluxはこの不定性の影響を大きく受けます。


3. XRTの望遠鏡の較正の不定性
XRTの有効面積レスポンス (arf) を計算するモンテカルロシミュレーション
 (xissimarfgen) 内では設計図どおりの理想的な望遠鏡の形状が仮定されてい
ます。しかし、実際には打ち上げ (+経年劣化?) によりXRTの形状がいくらか変化し
ていることが分かっており、XIS 0,1,3間、およびPINとのnormalizationの
食い違いの原因になります。この不定性はXIS 0,1,3で異なり、更に光軸(視野)
中心から離れるにつれて大きくなります。

4. 衛星の姿勢の揺れの不定性
観測中に「すざく」衛星の姿勢が揺れ、それが補正できないために天体が
光軸中心からずれてしまい、有効面積の見積もりを間違えてしまうことが
わかっています。XIS 0,1,3間、およびPINとのnormalizationの食い違いの
原因になります。この不定性は光軸(視野)中心から離れるにつれて大き
くなります。これまで割と良く補正されていたのですが、最近になってまた
姿勢の揺れが問題になりつつあるようです。aeattcor2というftoolを使うと少し
改善されるかもしれません。

2011年7月13日水曜日

FinderからFitsイメージをSnow Leopard用ds9 v6.2 X11で左クリックで開けるようにする

v6.1のaqua版を使え、というツッコミは正論すぎるので却下。
v6.1 aquaは何か不具合が生じたのでv6.2を使う。

1.  Automaterで「アプリケーション」を作成。
2.「シェルスクリプトを実行」をドラック&ドロップ
3.  以下みたいな感じ。-fileを付けないでうまく行かず詰まっていた。

4. 保存してApplicationフォルダに突っ込む。
(4+. できたAppを右クリックで「情報」を見て、一番上のアイコンのところにAqua版の”Appそのもの”をドラック&ドロップするとアイコンを変えることができる)
5. Fitsイメージを右クリックして関連付けを変更

2011年7月2日土曜日

Declination angleとEuler angle2の関係

いつもいつも忘れるが
EA2=90°-Dec

2011年6月15日水曜日

Astro-phに投稿した時にEPSファイルがうまく表示されない時

割とアホだが

1. Previewで開いてjpgで保存
2. jpgを再びPreviewで開いてPDF->PSで保存
3. ps2epsでEPSに変換

でうまくいった。

2011年5月28日土曜日

SciPyで数値積分 scipy.integrate.quad

#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);

2011年5月25日水曜日

Tabを4スペースに変換

アホみたいだが

sed 's/(TABキー出力)/(4スペース)/g' $1

と書いたスクリプトを使う。


追記: Eclipseなら「ソース→タブをスペースに変換」でできる。

2011年5月24日火曜日

BloggerでPythonソースのシンタックスハイライト

クリボウの Blogger Tips: コードをハイライトする「Code Prettify」ウィジェット
上記の方のウィジェットでできます。
<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()

Guplot.pyでデータの図をプロット

Gnuplot.pyをプロッタとして使用。
使い慣れている分、確かに使い易い。

#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')

出来た図は下。割と綺麗。



impotしたモジュールの元ファイル位置(インストールディレクトリ)を表示

__file__を使う。

In [515]: Gnuplot.__file__
Out[515]: '/Library/Frameworks/EPD64.framework/Versions/7.0/lib/python2.7/site-packages/Gnuplot/__init__.pyc'

のような感じ。

Pythonでモジュール一覧

自分世界 やってみるのが第一歩::Pythonでインストール済みモジュール一覧を表示

help('modules')

2011年5月22日日曜日

Pylabでデータの図をプロット

要はGnuplot+bashスクリプトからの卒業を目指して。

#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')


読んだデータファイルはこれ
出来た図は下。


しかし、pylabを使ってしまっているせいで結局どのオブジェクトについて扱っているのか、分かりにくくなってしまっている。
このサイトを参考に書きなおしてみる。
あと、ここをみると使いなれたGnuplotを使うのも悪くなさそう。
しかし、一番の目的であるFittingに到達する前に挫折。やれやれ。

Pylabの凡例 (legend)について

matplotlib artists — Matplotlib v1.0.1 documentation

defaultがnumpoints=2(1でない)なことにキレそうになった時にどうぞ。

2011年5月20日金曜日

ファイルからデータを読んでlistに突っ込む

for line in test_file:
    temp_list = float(line.strip().split(',')[0])
    test_list.append(temp_list)

Python で型を調べる

Python で型を調べる – types モジュール | すぐに忘れる脳みそのためのメモ

type(object)

でオブジェクトの型を返してくれる。
サンプルコード勉強中等は激しく便利。
良くも悪くも宣言しない時は何やっているのか訳が分からなくなる。

Python/インタラクティブシェルでの操作をファイルに保存する方法

Python/インタラクティブシェルでの操作をファイルに保存する方法 - TOBY SOFT wiki:
より

①前に入力した分も含めてまるごと保存

import readline
readline.write_history_file('test.py')

これはこれで便利

②IPythonのログ機能を使う方法
XSPCEライク?にデータを保存して行く方法

%logstartであらかじめログファイルを指定

%logstart 'test.py'
%logstopでログ記録終了

%logstop
%logon,%logoff でログのオン、オフ切り替え。


2011年4月12日火曜日

LEMOコネクタ? QLAコネクタ?

LEMOコネクタは商品名で所謂「LEMOコネクタ」以外のLEMOコネクタも存在するので正式名称を調べてみたのだが分からず。QLAコネクタも正式名称というよりは商品名らしいし…うーむ…。

2011年4月11日月曜日

GUIでネットリストを作成 on Mac OS X

GEDAというツールセットのgschemというコマンドを使う
GUIでネットリストを作成できるようです。



http://ryusai-hp.web.infoseek.co.jp/gaf.htm
下の方に「次のページ 回路図エディタgschem」とかがあるので
順番に進んでみて下さい。

finkでgedaをinstallするには

fink install geda-gaf

で簡単にできます。

2011年1月28日金曜日

XISRMFPARAM

サワDに対するメールより

ae_xi?_rmfparam_*.fits
の中を見るとParam_s*_p* という32行のパラメータ列
が並んでいると思います。これらがレスポンスの
詳細情報(メインピークの幅とメインピークに対
する他コンポーネントの幅・強度)を示しています。
各行の意味については
src/function/xisRespUtil/xisRespUtil.c
を読んで下さい。(これ以外ドキュメントがなさ
そうです・・・)

現在はメインピークの幅のみを変更しています。

Param_s*_p* の
2行目=s12
3行目=s13
4行目=s14
と定義してE(keV)に対して

FWHM(eV)=(S12+S13*E+S14*E**2)**0.5

です。(Param_s*_p* の1行目は1/(3.65*2.35)でPIとSigmaに変換するための係数
です)Eの単位はkeV、計算結果はeVなので注意して下さい。

なぜ、Param_s*_p*が複数あるかというと
s*は各セグメント毎にp*は同一セグメント内の場所(4カ所)毎に
異なったレスポンスを設定できるからです。

現在(ae_xi?_rmfparam_20080901.fits)は同一センサー
内では16列に全て同じ値を詰めています。

2011年1月12日水曜日

Xspecで異なるレスポンスかけてfit

diffuseや暗い点源の解析をやっているとCXBのモデルには異なるレスポンス(arf)をかけたくなることがあるが、Xspecでもそれは出来るようだ。
実際にやってみないと分からないが取りあえずマニュアルの対応箇所をペーストしておく。
変なマークアップが残ったままだけどまあいいか・・・

2011年1月6日木曜日

MITによるSpaced-row Charge Injection (SCI) の本論文


みんな Uchiyama et al. (2007) or (2009) を引いてくれるのを嬉しいのだが、あれはsaw-tooth correctionの論文であって、「すざく」の機上SCI論文としては下記を引くべき。


Prigozhin, Gregory; Burke, Barry; Bautz, Marshall; Kissel, Steve; Lamarr, Beverly
IEEE Transactions on Electron Devices, vol. 55, issue 8, pp. 2111-2120 (2008)

\bibitem[Prigozhin et al.(2008)]{2008ITED...55.2111P} Prigozhin, G., Burke, B., Bautz, M., Kissel, S.,  Lamarr, B.\ 2008, IEEE Transactions on Electron Devices, 55, 2111

MITの方が主著者で、査読付きのIEEE論文になっている、「筋としては」SCIの機上運用の本論文として引くべき一番の論文。…CIそのものについてが主でSCIの記述は少ないけど。


Bautz, M. W.; LaMarr, B. J.; Miller, E. D.; Kissel, S. E.; Prigozhin, G. Y.; Burke, B. E.; Gregory, J. A.; Uchiyama, H.; Hyodo, Y.; Yamaguchi, H.; Mori, H.; Tsuru, T.; Matsumoto, H.; Koyama, K.; Torii, K.; Katsuda, S.; Hasuike, K.; Nakajima, H.; Hayashida, K.; Tsunemi, H.; Murakami, H.; Ozaki, M.; Dotani, T.
UV, X-Ray, and Gamma-Ray Space Instrumentation for Astronomy XV. Edited by Siegmund, Oswald H. Proceedings of the SPIE, Volume 6686, pp. 66860Q-66860Q-11 (2007)

\bibitem[Bautz et al.(2007)]{2007SPIE.6686E..23B} Bautz, M.~W., LaMarr,
B.~J., Miller, E.~D., et al.\ 2007, \procspie, 6686

「すざく」でのSCI機上運用についての情報、見せたい(見せるべき)図などがまとまったBautz先生のSPIE論文。SCIの機上運用の論文としては内容的に最適。学生さんにも読んでもらいたいところ。…SPIE論文な上、Astro-phにも上がっていないので入手しにくいのが玉に瑕…。


でもUchiyama et al.も引いてくれると嬉しい。
\bibitem[Uchiyama et al.(2009)]{2009PASJ...61S...9U} Uchiyama, H., et al. 2009, \pasj, 61, 9