2017年12月11日月曜日

人間の耳で音の位相は聞き取れるのか?

「人間の耳で音の位相は聞き取れるのか?」問題のテスト。音痴な自分でも結構分かる気はする...(特にphase=0と0.5の差は相当顕著)。
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();
}

2017年12月7日木曜日

Pyorbitalで方位仰角、Rangeの計算

Pyorbitalで方位仰角、Rangeの計算。こんなのがさっくりできるのは、とてもありがたい...。どのくらい信用できるのかはわからないけど。
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