circle range

from astropy import units as u

from pylab import *

import numpy as np

import math as m

# test

def loadSatOrbitDat(datDir='',fileNum=50):

    oData = loadtxt(datDir+'1.txt')

    return oData

def transRaDec2D(ra, dec):

    x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795);

    y1 = np.cos(dec / 57.2957795) * np.sin(ra / 57.2957795);

    z1 = np.sin(dec / 57.2957795);

    return np.array([x1, y1, z1])

def locateSat(time=2459766., OrbitData = None, startId = 0, orbDataLen = 10000):

    pointfn.close()

def producePointingList3(out = '' , center = [60,-40], radius = 5, survey_file = 'skyMapOrSurveyList/E17.5_b17.5_beta_11.6_opt_transtime_1_CMG_1_dp_2_0.25_da_10_Texp_1.5_DEC60_500_0.1_800_1000_+5deg.dat'):

    points = np.loadtxt(survey_file)

    num = points.shape[0]

    center = center#ra dec

    radius = radius # degree

    c_eclip = coord.SkyCoord(points[:,2]*u.degree, points[:,1]*u.degree,frame='barycentrictrueecliptic')

    c_equtor = c_eclip.transform_to('icrs')

    cent_xyz = transRaDec2D(center[0],center[1])

    xyz = transRaDec2D((c_equtor.ra*u.degree).value,(c_equtor.dec*u.degree).value)

    xyz_dp = xyz[0]*cent_xyz[0] + xyz[1]*cent_xyz[1] + xyz[2]*cent_xyz[2]

    ids = xyz_dp>=m.cos(radius*m.pi/180.)

    # ids1 = (c_equtor.ra*u.degree).value > center[0]-radius

    # ids2 = (c_equtor[ids1].ra*u.degree).value < center[0]+radius

    # ids3 = (c_equtor[ids1][ids2].dec*u.degree).value > center[1]-radius

    # ids4 = (c_equtor[ids1][ids2][ids3].dec*u.degree).value < center[1]+radius

    obDataDir = 'orbit20160925/'

    orbitDat = loadSatOrbitDat(obDataDir,50)

    sOrbitId = 0

    elip_txt = '# ra    dec    lon(ecliptic)    lat(ecliptic)    pos_angle    time(julian)    sat_x    sat_y    sat_z    sun_x    sun_y    sun_z    moon_x    moon_y    moon_z    sat_vx    sat_vy    sat_vz    exp_time    isDeep\n'

    for p,p_ in zip(points[ids],c_equtor[ids]):

        ra = (p_.ra*u.degree).value

        dec = (p_.dec*u.degree).value

        # print(ra, dec)

        lon = p[2]

        lat = p[1]

        elip_txt = elip_txt + str(round(ra,6)) + '    ' + str(round(dec,6)) + '    ' + str(round(lon,6)) + '    ' + str(round(lat,6)) + '    -113.4333'

        oTime = p[0]

        satPos, sOrbitId =locateSat(time=oTime, OrbitData = orbitDat, startId = sOrbitId, orbDataLen = orbitDat.shape[0])

        tempOrbitId = sOrbitId

        deltT = 0.1 # unit s

        nTime = p[0]+deltT/86400.

        satPosN, _ =locateSat(time=nTime, OrbitData = orbitDat, startId = sOrbitId, orbDataLen = orbitDat.shape[0])

        sat_v = (satPosN - satPos)/deltT

        mm1 = np.sqrt(np.sum(sat_v*sat_v))

        mm2 = np.sqrt(np.sum(satPos*satPos))

        # print(np.sum(sat_v*satPos)/(mm1*mm2))

        # print(sat_v,p[6:9])

        if(np.abs(np.sum(sat_v*satPos)/(mm1*mm2))>0.01):

            print(oTime)

            break

        dat_col = [0,3,4,5,6,7,8,9,10,11,16,14]

        for col in dat_col[:-2]:

            elip_txt = elip_txt + '    ' + str(p[col])

        elip_txt = elip_txt + '    ' + str(sat_v[0]) + '    ' + str(sat_v[1]) + '    ' + str(sat_v[2])

        for col in dat_col[-2:]:

            elip_txt = elip_txt + '    ' + str(p[col])

        elip_txt = elip_txt + '\n'

    pointfn = open(out,'w')

    pointfn.write(elip_txt)

    pointfn.flush()

    pointfn.close()

if __name__ == "__main__":

    isRealSurvey = False

    isRealSurvey = True

    survey_file = 'skyMapOrSurveyList/E17.5_b17.5_beta_11.6_opt_transtime_1_CMG_1_dp_2_0.25_da_10_Texp_1.5_DEC60_500_0.1_800_1000_+5deg.dat'

    outFileName = 'pointing_test_false.dat'

    radi = 1.3

    outFileName = 'pointing_test_NGP_1..dat'

    radi = 1.2

    center_pos = [60, -40]

    center_pos = [192.8595, 27.1283]

    if isRealSurvey:

        producePointingList2(out = outFileName,center = center_pos, radius = radi, survey_file = survey_file)

        producePointingList3(out = outFileName,center = center_pos, radius = radi, survey_file = survey_file)

    else:

        producePointingList(out = outFileName,center = center_pos, radius = radi)