spacecruft
/
stvid
1
0
Fork 0
Code Releases Activity

Implemented acquisition on sunset to sunrise

pull/1/head
Cees Bassa 2018-03-10 14:09:16 +01:00
parent b47057ad82
commit 71940a9a4c
3 changed files with 144 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

49
acquire.py
View File

@ -1,13 +1,16 @@
#!/usr/bin/env python
from __future__ import print_function
import sys
import numpy as np
import cv2
import time
import ctypes
import multiprocessing
from astropy.coordinates import EarthLocation
from astropy.time import Time
from astropy.io import fits
from sun import get_sunriseset
import astropy.units as u
from utils import get_sunset_and_sunrise
# Capture images
def capture(buf,z1,t1,z2,t2,device,nx,ny,nz,tend):
@ -137,20 +140,37 @@ def compress(buf,z1,t1,z2,t2,nx,ny,nz,tend):
# Main function
if __name__ == '__main__':
# Current time
tnow=float(time.time())
tnow=1520573400.0
tnow=1520510226.281413,
tnow=Time.now()
# Set location
loc=EarthLocation(lat=52.8344*u.deg,lon=6.3785*u.deg,height=10*u.m)
# Reference altitude
refalt=-6.0*u.deg
# Get sunrise and sunset times
trise,tset=get_sunriseset(tnow,52.8344,6.3785,10,-6.0)
state,tset,trise=get_sunset_and_sunrise(tnow,loc,refalt)
print(tnow,trise,tset)
# Wait until sunset
if (tset<trise) & (tnow<tset):
print("Waiting for sunset at %s"%time.strftime("%FT%T",time.gmtime(tset)))
sys.exit()
# Start/end logic
if state=="sun never rises":
print("The sun never rises. Exiting program.")
sys.exit()
elif state=="sun never sets":
print("The sun never sets.")
tend=tnow+24*u.h
elif (trise<tset):
print("The sun is below the horizon.")
tend=trise
elif (trise>=tset):
dt=np.floor((tset-tnow).to(u.s).value)
print("The sun is above the horizon. Waiting %.0f seconds."%dt)
tend=trise
try:
time.sleep(dt)
except KeyboardInterrupt:
sys.exit()
print("Starting data acquisition. Acquisition will end at "+tend.isot)
# Settings
devid=int(sys.argv[1])
@ -176,12 +196,9 @@ if __name__ == '__main__':
t2=np.ctypeslib.as_array(t2base.get_obj())
buf=multiprocessing.Value('i',0)
tend=float(time.time())+31.0
# Set processes
pcapture=multiprocessing.Process(target=capture,args=(buf,z1,t1,z2,t2,device,nx,ny,nz,tend))
pcompress=multiprocessing.Process(target=compress,args=(buf,z1,t1,z2,t2,nx,ny,nz,tend))
pcapture=multiprocessing.Process(target=capture,args=(buf,z1,t1,z2,t2,device,nx,ny,nz,tend.unix))
pcompress=multiprocessing.Process(target=compress,args=(buf,z1,t1,z2,t2,nx,ny,nz,tend.unix))
# Start
pcapture.start()

77
plot.py
View File

@ -1,69 +1,30 @@
#!/usr/bin/env python
import sys,os,glob
from stio import fourframe,satid,observation
import numpy as np
import sys
from astropy.io import fits
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fname="2018-02-12T17:39:55.270.fits"
ff=fourframe("2018-02-26T03:09:05.866.fits")
# Read fits
hdu=fits.open(fname)
# Read data
data=hdu[0].data
zavg,zstd,zmax,znum=data
# Read header
ny,nx=zavg.shape
nz=hdu[0].header['NFRAMES']
dt=[hdu[0].header['DT%04d'%i] for i in range(nz)]
# Sigma frame
zsig=(zmax-zavg)/(zstd+1e-9)
# Position
xm,ym=np.meshgrid(np.arange(nx),np.arange(ny))
# Selection width
w=10
zsig=(ff.zmax-ff.zavg)/(ff.zstd+1e-9)
sigma=5.0
# Selection
c=(zsig>sigma) & (xm>w) & (xm<nx-w) & (ym>w) & (ym<ny-w)
x=np.ravel(xm[c])
y=np.ravel(ym[c])
inum=np.ravel(znum[c]).astype('int')
c=zsig>sigma
print(np.sum(c)/float(ff.nx*ff.ny))
xm,ym=np.meshgrid(np.arange(ff.nx),np.arange(ff.ny))
x,y=np.ravel(xm[c]),np.ravel(ym[c])
inum=np.ravel(ff.znum[c]).astype('int')
sig=np.ravel(zsig[c])
t=np.array([dt[i] for i in inum])
t=np.array([ff.dt[i] for i in inum])
# Load predictions
d=np.loadtxt(fname+".id",usecols=(1,2,3,4,5,6))
plt.figure()
plt.scatter(x,y,c=t,s=1)
plt.colorbar()
plt.show()
x0=d[:,0]
y0=d[:,1]
t0=np.zeros_like(x0)
x1=d[:,2]
y1=d[:,3]
t1=d[:,4]
for i in range(len(x0)):
xr=x0[i]+(x1[i]-x0[i])*(t-t0[i])/(t1[i]-t0[i])
yr=y0[i]+(y1[i]-y0[i])*(t-t0[i])/(t1[i]-t0[i])
r=np.sqrt((x-xr)**2+(y-yr)**2)
c=r<10.0
print(i,np.sum(c))
plt.figure()
plt.plot([x0[i],x1[i]],[y0[i],y1[i]],'r')
plt.scatter(x[c],y[c],c=t[c],s=10.0*(sig[c]-sigma))
plt.scatter(x,y,c=t,s=1)
plt.show()
#fig=plt.figure()
#ax=fig.add_subplot(111,projection='3d')
#ax.scatter(x,y,t,c=t,s=10*(sig-5))
#for i in range(len(x0)):
# ax.plot([x0[i],x1[i]],[y0[i],y1[i]],[t0[i],t1[i]],'k')
#plt.show()
fig=plt.figure()
ax=fig.add_subplot(111,projection='3d')
ax.scatter(x,y,t,c=t,s=1)
plt.show()

92
utils.py 100644
View File

@ -0,0 +1,92 @@
from astropy.coordinates import SkyCoord,EarthLocation,AltAz,FK5,get_sun
from astropy.time import Time
import astropy.units as u
import numpy as np
from scipy import interpolate
# Sunrise/sunset algorithm from Astronomical Algorithms by Jean Meeus
def get_sunset_and_sunrise(tnow,loc,refalt):
# Get time
nmjd=64
mjd0=np.floor(tnow.mjd)
mnow=tnow.mjd-mjd0
mjd=np.linspace(mjd0-1.0,mjd0+2.0,nmjd)
t=Time(mjd,format='mjd',scale='utc')
# Get sun position
psun=get_sun(t)
# Correct for precession by converting to FK5
pos=SkyCoord(ra=psun.ra.degree,dec=psun.dec.degree,frame='icrs',unit='deg').transform_to(FK5(equinox=t))
# Compute altitude extrema
de=np.mean(pos.dec)
minalt=np.arcsin(np.sin(loc.latitude)*np.sin(de)-np.cos(loc.latitude)*np.cos(de))
maxalt=np.arcsin(np.sin(loc.latitude)*np.sin(de)+np.cos(loc.latitude)*np.cos(de))
# Never sets, never rises?
if minalt>refalt:
return "sun never sets",t[0],t[0]
elif maxalt<refalt:
return "sun never rises",t[0],t[0]
# Set up interpolating function
fra=interpolate.interp1d(t.mjd,pos.ra.degree)
fde=interpolate.interp1d(t.mjd,pos.dec.degree)
# Get GMST
gmst0=Time(mjd0,format='mjd',scale='utc').sidereal_time('mean','greenwich')
# Get transit time
mtransit=np.mod((fra(mjd0)*u.degree-loc.longitude-gmst0)/(360.0*u.degree),1.0)
while True:
gmst=gmst0+360.985647*u.deg*mtransit
ra=fra(mjd0+mtransit)*u.deg
ha=gmst+loc.longitude-ra
mtransit-=ha/(360.0*u.deg)
if np.abs(ha.degree)<1e-9:
break
# Set transit
ttransit=Time(mjd0+mtransit.value,format='mjd',scale='utc')
# Hour angle offset
ha0=np.arccos((np.sin(refalt)-np.sin(loc.latitude)*np.sin(np.mean(pos.dec)))/(np.cos(loc.latitude)*np.cos(np.mean(pos.dec))))
# Get rise time
mrise=mtransit-ha0/(360.0*u.deg)
if mrise<mnow:
mrise+=1.0
while True:
gmst=gmst0+360.985647*u.deg*mrise
ra=fra(mjd0+mrise)*u.deg
de=fde(mjd0+mrise)*u.deg
ha=gmst+loc.longitude-ra
alt=np.arcsin(np.sin(loc.latitude)*np.sin(de)+np.cos(loc.latitude)*np.cos(de)*np.cos(ha))
dm=(alt-refalt)/(360.0*u.deg*np.cos(de)*np.cos(loc.latitude)*np.sin(ha))
mrise+=dm
if np.abs(dm)<1e-9:
break
# Set rise time
trise=Time(mjd0+mrise.value,format='mjd',scale='utc')
# Get set time
mset=mtransit+ha0/(360.0*u.deg)
if mset<mnow:
mset+=1.0
while True:
gmst=gmst0+360.985647*u.deg*mset
ra=fra(mjd0+mset)*u.deg
de=fde(mjd0+mset)*u.deg
ha=gmst+loc.longitude-ra
alt=np.arcsin(np.sin(loc.latitude)*np.sin(de)+np.cos(loc.latitude)*np.cos(de)*np.cos(ha))
dm=(alt-refalt)/(360.0*u.deg*np.cos(de)*np.cos(loc.latitude)*np.sin(ha))
mset+=dm
if np.abs(dm)<1e-9:
break
# Set set time
tset=Time(mjd0+mset.value,format='mjd',scale='utc')
return "sun rises and sets",tset,trise