Aide-memoire: 2016

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

登録: 投稿 (Atom)

Aide-memoire

2016年11月26日土曜日

Processingでブラウン運動

2016年10月5日水曜日

2HD with Arduno

2016年10月3日月曜日

2HD

2016年5月23日月曜日

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

2016年5月10日火曜日

Labviewでジャンケン

2016年4月19日火曜日

お察し下さい

2016年4月7日木曜日

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

2016年3月31日木曜日

Arf_BG_corrector.py

2016年3月27日日曜日

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

2016年2月24日水曜日

PI Shifter for Suzaku XIS

プロフィール

ブログアーカイブ

ラベル

このブログを検索

2016年11月26日土曜日

2016年10月5日水曜日

2016年10月3日月曜日

2016年5月23日月曜日

2016年5月10日火曜日

2016年4月19日火曜日

2016年4月7日木曜日

2016年3月31日木曜日

2016年3月27日日曜日

2016年2月24日水曜日

プロフィール

ブログ アーカイブ

ラベル

このブログを検索

ブログアーカイブ