/ This function calculate XYZ coordinate of a GPS satellite in a RINEX file

/ Sat clock error(NANO SECOND) stored as SD value of the return IdPt3D object
/ Convert to second by multiplying with 1E-9

/ Author : anonymous

/ ***************************************************************************

  Function IdPt3D XYZsat = rinex2xyz(String fname, Date t_date, int prn_no)

/ ***************************************************************************

  dow = t_date.dow()
  tm  = t_date.time()
  t   = (dow-1)*24*3600 + tm*3600

  V = VecStr()
  V.load(fname)

  n = V.size()
  line_no = 0

/ Start searching for the requiring PRN NO at the 4th line
  for i=3, n-1
    no = int(V(i).left(2).val())
    if no==prn_no
      line_no = i
      break
    end
  end

if line_no==0
  print "PRN no not found"
  return
end

if line_no+7 > n-1
  print "ERROR : no of data lines are not enough for the satellite"
  return
end


/ Read line 1

/ ***************** PRN / EPOCH / SV CLK *******************

line   = V(line_no).left(79)

prn    = int(line.sub(0,1).val())
year   = int(line.sub(2,4).val())
month  = int(line.sub(5,7).val())
day    = int(line.sub(8,10).val())
hour   = int(line.sub(11,13).val())
minute = int(line.sub(14,16).val())
sec    = line.sub(17,21).val()

clk_bias   = line.sub(22,40).val()
clk_drift  = line.sub(41,59).val()
clk_driftr = line.sub(60,78).val()

/ ***************** BROADCAST ORBIT - 1 *******************

line    = V(line_no + 1).left(79)

IODE    = line.sub(3,21).val()
Crs     = line.sub(22,40).val()
Delta_n = line.sub(41,59).val()
M0      = line.sub(60,78).val()

/ ***************** BROADCAST ORBIT - 2 *******************

line     = V(line_no + 2).left(79)

Cuc      = line.sub(3,21).val()
e        = line.sub(22,40).val()
Cus      = line.sub(41,59).val()
square_A = line.sub(60,78).val()

/ ***************** BROADCAST ORBIT - 3 *******************

line   = V(line_no + 3).left(79)

toe    = line.sub(3,21).val()
Cic    = line.sub(22,40).val()
OMEGA0 = line.sub(41,59).val()
Cis    = line.sub(60,78).val()

/ ***************** BROADCAST ORBIT - 4 *******************

line      = V(line_no + 4).left(79)

i0        = line.sub(3,21).val()
Crc       = line.sub(22,40).val()
omega     = line.sub(41,59).val()
OMEGA_DOT = line.sub(60,78).val()

/ ***************** BROADCAST ORBIT - 5 *******************

line     = V(line_no + 5).left(79)

IDOT     = line.sub(3,21).val()
L2_code  = line.sub(22,40).val()
week_no  = line.sub(41,59).val()
L2_pflag = line.sub(60,78).val()

/ ***************** BROADCAST ORBIT - 6 *******************

line      = V(line_no + 6).left(79)

SV_acc    = line.sub(3,21).val()
SV_health = line.sub(22,40).val()
TGD       = line.sub(41,59).val()
IODC      = line.sub(60,78).val()

/ ***************** BROADCAST ORBIT - 7 *******************

line   = V(line_no + 7).left(79)

HOW    = line.sub(3,21).val()
spare1 = line.sub(22,40).val()
spare2 = line.sub(41,59).val()
spare3 = line.sub(60,78).val()


/ This program compute Position and velocity of satellite from
/ 6 Kleperian elements
set angle "rad"
set format "fix"
set precision 10


/a    = orbit semi-major axis
/e    = orbit eccentricity
/i    = orbit inclination angle
/rasc = right ascension of node
/argm = argument of perigee
/M    = mean anomaly
/we   = WGS84 value of the earth rotation rate

we     = 7.2921151467e-5
GM     = 3.986005E14
TRESH  = 0.00000000000001
MAX_no = 30

tk = t - toe

a  = square_A ^2

n0 = sqrt(GM/a^3)
n  = n0 + Delta_n
M  = M0+n*tk

/ Calculate Eccentric anomaly

E1 = M

diff = 100

no = 0
while (diff > TRESH & no<MAX_no)
   E2   = M + e * sin(E1)
   diff = abs(E2 - E1)
   E1   = E2
   no   = no+1
end

Ec = E2

/ Calculate True anomaly (nu)

temp = (1-e^2)^0.5

nu = atan2(temp * sin(Ec), (cos(Ec) - e))

PHE = nu + omega

Cor_uk = Cuc*cos(2*PHE) + Cus*sin(2*PHE)
Cor_rk = Crc*cos(2*PHE) + Crs*sin(2*PHE)
Cor_ik = Cic*cos(2*PHE) + Cis*sin(2*PHE)

uk = PHE + Cor_uk
rk = a*(1-e*cos(Ec)) + Cor_rk
ik = i0 + IDOT*tk + Cor_ik

Xpk = rk*cos(uk)
Ypk = rk*sin(uk)

OMEGA = OMEGA0 + (OMEGA_DOT - we)*tk - we*toe

Xk = Xpk*cos(OMEGA) - Ypk*sin(OMEGA)*cos(ik)
Yk = Xpk*sin(OMEGA) + Ypk*cos(OMEGA)*cos(ik)
Zk = Ypk*sin(ik)

/ ****** The following is to calculate satellite clock error
/ unit in NANO SEC

if year > 80
   year = year + 1900
else
   year = year + 2000
end

date0 = Date(year, month, day, hour+minute/60+sec/3600)

jd0 = date0.jd()
jd1 = t_date.jd()
dt  = (jd1-jd0)*24*3600
dt_error = (clk_bias + clk_drift*dt + clk_driftr*dt^2)*1e9

XYZsat   = IdPt3D(prn_no, Xk, Yk, Zk, dt_error)

return
 

/ This function returns a matrix of size n by n, of a cross symbol
/ The cross is one pixel thick and will have a value = value2
/ The background of the matrix has a value = value1

/ Author ; anonymous

/ *********************************************************************

  Function Matrix out = sym_cross(int size, double value1, double value2)

/ *********************************************************************

/ value1 is color of background
/ value2 is color of the cross

  if is_even(size)
    print "Error in external function sym_cros : size not an odd number"
    return
  end

  out = Matrix(size, size, value1)
  mid = fix(size/2)

  for i=0,size-1
    out(mid, i) = value2
    out(i, mid) = value2
  end

  return