Minor Planet Center files to Celestia

[amaury]

Please do not use the current script as its output for comets is inaccurate. However the asteroid part is ok. I will upload a better one when it is updated.
Amaury


Original post below:
-----------------------------------------------------
my scripting frenzy is not getting any better

i'm making a celx script to convert .dat files from the Minor Planet Center to scc files :

http://asolignac.free.fr/celestia/MPC_to_Celestia_translator.celx

it requires
ftp://cfa-ftp.harvard.edu/pub/MPCORB/NEAcr.dat
and/or
ftp://cfa-ftp.harvard.edu/pub/MPCORB/COMET.dat
(anonymous ftp login, or you can try http://cfa-www.harvard.edu/iau/MPCORB.html for mirrors)

the resulting files are an update of Arlene and Dirl's motherlode files, EXCEPT i'm still struggling with the radius from magnitude estimation:
I get giant comets!!

any hint, Selden ?

Amaury

[selden]

Well, I get a "no such file" error

You aren't trying to apply the algorithm used for asteroids to a comet, are you? Asteroids are much darker.

[symaski62]

my scripting frenzy is not getting any better

i'm making a celx script to convert .dat files from the Minor Planet Center to scc files :

http://asolignac.free.fr/celestia/MPC_to_Celestia_translator.celx

it requires
ftp://cfa-ftp.harvard.edu/pub/MPCORB/NEAcr.dat
and/or
ftp://cfa-ftp.harvard.edu/pub/MPCORB/COMET.dat
(anonymous ftp login, or you can try http://cfa-www.harvard.edu/iau/MPCORB.html for mirrors)

the resulting files are an update of Arlene and Dirl's motherlode files, EXCEPT i'm still struggling with the radius from magnitude estimation:
I get giant comets!!

any hint, Selden ?

Amaury

it requires
ftp://cfa-ftp.harvard.edu/pub/MPCORB/NEAcr.dat
and/or
ftp://cfa-ftp.harvard.edu/pub/MPCORB/COMET.dat

no url



charge

http://mpcorb.astro.cz/
COMET.DAT
DAILY.DAT



\MPC\MPC_to_Celestia_translator.celx
\MPC\COMET.dat
\MPC\DAILY.dat

MPC_to_Celestia_translator.celx "CLIC"

Escape => ERROR :/

MPC_to_Celestia_translator_V1.1.celx

[amaury]

\MPC\MPC_to_Celestia_translator.celx
\MPC\COMET.dat
\MPC\DAILY.dat

MPC_to_Celestia_translator.celx "CLIC"

Escape => ERROR :/

MPC_to_Celestia_translator_V1.1.celx

Uh ?!?
do you get an error when you press escape ?

[amaury]

i get it : there should be a wait() right after c:requestsystemaccess() on line 242...

i didn't get the error myself because my celestia config file had systemaccess on auto (i got tired of pressing (y) all the time while testing )

i can't upload anything from work, and no home connection for a week at least, so you'll have to insert yourself the wait() command between line 242 and 243 if you want to test again.

If anyone wants to modify & upload in this thread (not on motherlode), feel free to.

amaury

thanks for testing, Symaski and Selden!

Selden : yes, you're right, i should find a proper comet albedo, and a more specific radius formula...

[amaury]

Still working on it... (wrong radius formula)

amaury

[Dirl]

Is there someone have this celx script? Amaury's link do not work, and he did not visit the forum since 2007.

[amaury]

Hi Dirl,

i haven't been cel-scripting for quite a while now, and the links for the scripts are down because the hosting company closed my account without notice last year.

I will look into my old computer files and post again if i find it

Amaury

[amaury]

Hello,

5 years later I found my stuff back (ouch!)

These are scripts adapted from the MPC Perl script by Arlene Ducao.

I made them a while ago, so do NOT launch them without reading the content and adjusting URLs, etc. Otherwise you are likely to get errors all over

Awful code by the way, but hey, we all learn on the way!

NB : do not trust comet albedo, it is not realistic.

Minor Planet Center >> Celestia

Code: Select all

--********************************************************************
--********************************************************************
--****** Minor Planet Center files to Celestia .SSC translator *******
--********************************************************************
--********************************************************************

-- Script adapted to Celx by Amaury Solignac, february 2006.
-- (original Perl script by Arlene Ducao)

-- Useful information gathered from Selden Ball's page on:
-- "How to transform orbital elements into Celestia's SSC format"
-- http://www.lns.cornell.edu/~seb/celestia/transforming_ephemeris.html
-- *******************************************************************
--
-- This script will turn .DAT files from the Minor Planet Center into
-- Celestia .SSC files.
--
-- It requires one or both of these files:
-- http://www.astro.cz/mpcorb/COMET.DAT
-- http://www.astro.cz/mpcorb/NEA.DAT
--
-- You may also download them from the Minor Planet Center FTP:
-- ftp://cfa-ftp.harvard.edu/pub/MPCORB/
--
-- The Period formula is different for hyperbolic vs. elliptic comets.
-- Comets with parabolic orbits will be turned into hyperbolic ones
-- by adding 0.00001 to their original eccentricity (1).
-- See http://shatters.net/forum/viewtopic.php?t=8940 for details.
--
-- While some values are calculated, others, such as albedo and radius,
-- can only be estimated.
-- Therefore:
--   Near Earth Asteroids are given an mean albedo of 0.11 which is used
--   to estimate the radius.
--   Comets are given an albedo of 0.03, and their nuclei radius is
--   estimated to a *very* rough average of 2km.
--
-- NB:
-- If you are running Windows and plan on using the .DAT files for
-- something else, you should also download NEAcr.dat and COMETcr.dat.
--
-- Nothing prevents you from modifying the script to run it with
-- the full MPCORB.dat file, but it produces a *HUGE* SSC file that's
-- probably not yet suited for Celestia on our home computers...
-- You will also have to remove a few comment lines at the beginning of
-- the file. Finally, the average NEA albedo value of 0.11 will not fit.
--
--*********************************************************************


-------------------------------------------------
-- FUNCTIONS ------------------------------------
-------------------------------------------------

function celestia_keyboard_callback(input)
------------------------------------------
   key_pressed = input
   return(true)
end
---


function convert(input_name, input_extension)
--------------------------------------------
   date_stamp = os.date("*t")
   date_stamp = string.format("%04d_%02d_%02d", date_stamp.year, date_stamp.month, date_stamp.day)

   num = 0
   output_file = input_name.."_MPC_"..date_stamp..".ssc"
   input_file = input_name.."."..input_extension

 if (io.open(input_file,"r")) then
    io.close()
   output = io.open(output_file,"w")
   output:write("# Data extracted from "..input_file.." (Minor Planet Center)\n\n\n")

    if input_name == "COMET" then
      for line in io.lines(input_file) do
         COMET_print(line, input_name)
         c:print(header..message.." Processing "..input_file..": "..num, 10, -1, 0)
         wait()
      end
    else
      for line in io.lines(input_file) do
         ASTEROID_print(line, input_name)
         c:print(header..message.." Processing "..input_file..": "..num, 10, -1, 0)
         wait()
      end
    end

   message = message.." "..num.. " objects written to "..output_file.."\n"
   output = io.close()
 else
    message = message.." ****Oops! Couldn't find "..input_file..": "..output_file.." not created...\n"

 end
end
---


function ASTEROID_print(line, input_name)
-----------------------------------------

------------- Parameter gathering

   H = string.sub(line, 9, 13)                           -- Absolute magnitude
   Epoch = string.sub(line, 21, 25)                     -- Epoch (in packed form, .0 TT)
                                                --[[
                                                The first two digits of the year are packed into a single character in column
                                                1 (I = 18, J = 19, K = 20). Columns 2-3 contain the last two digits of the year.
                                                Column 4 contains the month and column 5 contains the day, coded as detailed
                                                below:
                                                   Month     Day      Character         Day      Character
                                                                in Col 4 or 5              in Col 4 or 5
                                                   Jan.       1           1             17           H
                                                   Feb.       2           2             18           I
                                                   Mar.       3           3             19           J
                                                   Apr.       4           4             20           K
                                                   May        5           5             21           L
                                                   June       6           6             22           M
                                                   July       7           7             23           N
                                                   Aug.       8           8             24           O
                                                   Sept.      9           9             25           P
                                                   Oct.      10           A             26           Q
                                                   Nov.      11           B             27           R
                                                   Dec.      12           C             28           S
                                                          13           D             29           T
                                                          14           E             30           U
                                                          15           F             31           V
                                                          16           G           
                                                
                                                Examples:
                                                   1996 Jan. 1    = J9611
                                                   1996 Jan. 10   = J961A
                                                   1996 Sept.30   = J969U
                                                   1996 Oct. 1    = J96A1
                                                   2001 Oct. 22   = K01AM
                                                ]]

      -- Let's unpack the epoch using the MPC packed format rules (see long comment between [[brackets]] above)
      if (string.sub(Epoch, 1, 1) == "I") then
         EpochC = "18"
      elseif (string.sub(Epoch, 1, 1) == "J") then
         EpochC = "19"
      else
         EpochC = "20"
      end
   EpochY = EpochC..string.sub(Epoch, 2, 3)
   EpochM = tonumber(string.sub(Epoch, 4, 4), 31)            -- uses a convenient LUA tonumber(arg, base) native function, with base > 10, to fit the rule below
   EpochD = tonumber(string.sub(Epoch, 5, 5), 31)
   Julian_Epoch = c:tojulianday(EpochY, EpochM, EpochD, 0, 0, 0)


   M = string.sub(line, 27, 35)                        -- Mean anomaly at the epoch, in degrees
   Peri = string.sub(line, 38, 46)                        -- Argument of perihelion, J2000.0 (degrees)
   Node = string.sub(line, 49, 57)                        -- Longitude of the ascending node, J2000.0 (degrees)
   Incl = string.sub(line, 60, 68)                        -- Inclination to the ecliptic, J2000.0 (degrees)
   e = string.sub(line, 71, 79)                        -- Orbital eccentricity
   a =  string.sub(line, 93, 103)                        -- Semimajor axis (AU)
   U = string.sub(line, 106, 106)                        --[[
                                                 Uncertainty parameter, U
                                                 If this column contains `E' it indicates
                                                that the orbital eccentricity was assumed.
                                                For one-opposition orbits this column can
                                                also contain `D' if a double (or multiple)
                                                designation is involved or `F' if an e-assumed
                                                double (or multiple) designation is involved.
                                                ]]


   Ref = string.sub(line, 108, 116)                     -- Reference

   orbit_type = string.sub(line, 162, 165)                  -- Not yet implemented in my script...
                                                --[[
                                                 (4-hexdigit flags)
                                                The bottom 6 bits are used to encode a
                                                value representing the orbit type (other
                                                values are undefined):
                                                  Value
                                                   2  Aten
                                                   3  Apollo
                                                   4  Amor
                                                   5  Object with q < 1.381 AU
                                                   6  Object with q < 1.523 AU
                                                   7  Object with q < 1.665 AU
                                                Additional information is conveyed by
                                                adding in the following bit values:
                                                16384  Critical list numbered object
                                                32768  Object is PHA
                                                ]]



   Designation = string.gsub(string.sub(line, 175, 194), " ", "")   -- Readable designation


   LastObs = string.sub(line, 195, 202)                  -- Date of last observation included in orbit solution (YYYYMMDD format)

   albedo = 0.11                                    -- a commonly assumed average value for NEAs

   radius = 0.5 * 1329 *((math.pow(10, -(H)/5))/math.sqrt(albedo))
      -- http://www.physics.sfasu.edu/astro/asteroids/sizemagnitude.html
      -- http://cfa-www.harvard.edu/iau/lists/Sizes.html
      
------------- Output      

      output:write("# Asteroid "..Designation.." from the Minor Planet Center "..input_file.." file\n") -- Unaltered designation

      Designation = string.gsub(Designation, " ", "")       -- To avoid spaces in the object's name
      Designation = string.gsub(Designation, "/", " ")       -- To avoid slashes in the object's name

      output:write("\""..Designation.."\" \"Sol\"\n")
      output:write("{\n")
      output:write("   Class \"asteroid\"\n")
      output:write("   Mesh \"asteroid.cms\"\n")
      output:write("   Texture \"asteroid.jpg\"\n")
      output:write("    Radius  "..radius.."          #  estimated: 0.5 * 1329 *(10^(-AbsMag[=H]/5))/sqrt(albedo))\n")
      output:write("   EllipticalOrbit\n   {\n")
      output:write("      Epoch  "..Julian_Epoch.."\n")
      output:write("      Period "..math.sqrt(math.pow(a,3)).."    #  calculated: sqrt(SemiMajorAxis**3)\n")
      output:write("      SemiMajorAxis "..a.." \n")
      output:write("      Eccentricity "..e.." \n")
      output:write("      Inclination "..Incl.." \n")
      output:write("      AscendingNode "..Node.." \n")
      output:write("      ArgOfPericenter "..Peri.." \n")
      output:write("      MeanAnomaly "..M.." \n")
      output:write("   }\n")
      output:write("   Albedo            "..albedo.."  # an average albedo value for Near Earth Asteroids\n")
      output:write("   RotationPeriod    "..math.random(365).."   # random value (0<RP<365)\n")
      output:write("}\n\n\n")

      num = num+1

end
---




function COMET_print(line, input_name)
--------------------------------------

--[[ From Selden Ball's page on how to transform orbital elements to Celestia's SSC formaat
   (http://www.lns.cornell.edu/~seb/celestia/transforming_ephemeris.html):
   
   "By convention, cometary elements don't include a mean anomaly value, so the Epoch you put into
   the SSC's EllipticalOrbit must be the time of the pericenter passage, when the mean anomaly is,
   by definition, 0.0."]]
   
------------- Parameter gathering

   M = 0                                          -- Mean anomaly at the time of perihelion, in degrees
   H = string.sub(line, 92, 95)                        -- Absolute magnitude
   q = string.sub(line, 31, 39)                        -- Perihelion distance
   PeriY = string.sub(line, 15, 18)                     -- Perihelion year
   PeriM = string.sub(line, 20, 21)                     -- Perihelion month
   PeriD = string.sub(line, 23, 24)                     -- Perihelion day
   Julian_PeriTime = c:tojulianday(PeriY, PeriM, PeriD, 0, 0, 0) -- Perihelion date in julian format for Celestia
   Peri = string.sub(line, 52, 59)                        -- Argument of perihelion, J2000.0 (degrees)
   Node = string.sub(line, 62, 69)                        -- Longitude of the ascending node, J2000.0 (degrees)
   Incl = string.sub(line, 72, 79)                        -- Inclination to the ecliptic, J2000.0 (degrees)
   e = tonumber(string.sub(line, 42, 49))                  -- Orbital eccentricity
   Designation = string.sub(line, 103, -1)                  -- Readable designation
   Designation = string.sub(line, 103, -1)                  -- Readable designation

   albedo = 0.03                                    -- Estimate of average albedo for active comets NUCLEI
      -- From "Extreme albedo comets and impact hazard" (2004) W. M. Napier, J. T. Wickramasinghe and N. C. Wickramasinghe

   radius = 2
      -- a *very* rough estimated average for comet nuclei radius
      -- the problem is that the nucleus radius cannot be directly related to the absolute magnitude of the comet, because of the chemical processes involved

   if e < 1 then
     period = (math.pow(q/(1-e), 1.5))                     -- for elliptic orbits
     period_message = "    # calculated: (q/(1-e))**1.5)\n"
    else
       if e == 1 then e = 1 + 0.00001 end                  -- to turn parabolic orbits into hyperbolic ones (to avoid div by 0)
     period = (math.pow(q/(e-1), 1.5))                     -- for parabolic / hyperbolic orbits
     period_message = "    # calculated: (q/(e-1))**1.5) (hyperbolic orbit)\n"
   end
   
------------- Output      

      output:write("# Comet "..Designation.." from the Minor Planet Center "..input_file.." file\n") -- Unaltered designation

      Designation = string.gsub(Designation, "/", "")       -- To avoid slashes in the object's name
      Designation = string.gsub(Designation, " ", "")       -- To avoid spaces in the object's name

      output:write("\""..Designation.."\" \"Sol\"\n")
      output:write("{\n")
      output:write("   Class \"comet\"\n")
      output:write("   Mesh \"halley.cmod\"\n")
      output:write("   Texture \"asteroid.jpg\"\n")
      output:write("    Radius  "..radius.."          #  a *very* rough average\n")
      output:write("   EllipticalOrbit\n   {\n")
      output:write("      PericenterDistance "..q.."\n")
      output:write("      Period "..period..period_message)
      output:write("      Eccentricity "..e.." \n")
      output:write("      Inclination "..Incl.." \n")
      output:write("      AscendingNode "..Node.." \n")
      output:write("      ArgOfPericenter "..Peri.." \n")
      output:write("      MeanAnomaly "..M.." \n")
      output:write("      Epoch  "..Julian_PeriTime.."\n")
      output:write("   }\n")
      output:write("   Albedo            "..albedo.."  # estimate of average albedo for active comets NUCLEI\n")   
      output:write("   RotationPeriod    "..math.random(365).."   # random value (0<RP<365)\n")
      output:write("}\n\n\n")
      num = num+1
end
---


-------------------------------------------------
-- SCRIPT START----------------------------------
-------------------------------------------------
script_version = "MPC .dat to Celestia .ssc converter (beta2)"
header = " "..script_version.."\n ______________________________________\n"
c = celestia


c:requestkeyboard(true)
while(key_pressed ~= " ") do
   c:print(header.." This script will look for .DAT files in your\n Celestia root folder and convert them to\n Celestia .SSC files.\n\n You will need at least one of these two files\n from the Minor Planet Center:\n http://www.astro.cz/mpcorb/COMET.DAT\n http://www.astro.cz/mpcorb/NEA.DAT\n\n If you already downloaded one or both files,\n press the Spacebar to run the conversion.\n Otherwise, press the Escape key to exit now...", 1, -1, 0)
   wait()
end
c:requestkeyboard(false)
c:print(header.." Celestia will now request system access\n in order to handle the output file(s).", 5, -1, 0)
wait(4)
c:requestsystemaccess()
wait()
num = 0
message = ""
input_extension = ""


convert("COMET", "dat")
convert("NEA", "dat")

c:print(header..message.."\n End of script.", 10, -1, 0)
wait(10)




--[[

Typical .ssc asteroid item :
----------------------------
"Ceres" "Sol"
{
   Class "asteroid"
   Texture "asteroid.jpg"
   Radius 487.5        # 975x909km
   Oblateness 0.068    # Thomas et al, Nature, 8 Sep 2005
   EllipticalOrbit
   {
      Epoch       2452600.5     # 2002 Nov 22 00:00UT
      Period            4.60014
      SemiMajorAxis     2.7660
      Eccentricity      0.0793
      Inclination      10.584
      AscendingNode    80.483
      ArgOfPericenter  74.043
      MeanAnomaly     232.067
   }
   RotationPeriod         9.075
   RotationEpoch    2449249.91125 # 1993 Sep 19 09:52:12UT
   Obliquity             11
   EquatorAscendingNode  29
   RotationOffset       339.85
   Albedo 0.113
}



Typical .ssc comet item :
----------------------------
"C1995 O1 (Hale-Bopp)" "Sol"
{
   Class "comet"
   Mesh "asteroid.cms"
   Texture "asteroid.jpg"
   Radius 50.0      # guess
        Albedo 0.1
   EllipticalOrbit
   {
      Period             2550.259099
      PericenterDistance    0.931444
      Eccentricity          0.995010
      Inclination          89.36570
      AscendingNode       282.50470
      ArgOfPericenter     130.85030
        MeanAnomaly           0.0
      Epoch           2450537.8292
   }
   InfoURL "http://ssd.jpl.nasa.gov/cgi-bin/eph"
}

]]


JPL Small Body Database >> Celestia

Code: Select all

--********************************************************************
--***** JPL Small Body Database files to Celestia .SSC translator ****
--                   Amaury Solignac (april 2006)

--*************************** version 1 ******************************
-- CREDITS:
-- Adapted from an original MPC Perl script by Arlene Ducao.
-- Useful information gathered from Selden Ball's page on:
-- "How to transform orbital elements into Celestia's SSC format"
-- http://www.lns.cornell.edu/~seb/celestia/transforming_ephemeris.html
-- Finally, the fromCSV function was adapted from Roberto Ierusalimschy
-- *******************************************************************


-- This script will turn .CSV files from the JPL SBDB into .SSC files
-- for use in celestia.
-- Visit http://ssd.jpl.nasa.gov/sbdb_query.cgi for details

-- The reasons for using the JPL SBDB rather than the JPL Horizons or
-- Minor Planet center files are that SBDB can contain both orbital
-- and physical data, and that elements accuracy is better.



-- ONCE YOU ARE READY to USE THE SCRIPT:

-- 1)
--  Copy and paste one of these URLs into your favorite browser:

--    ## for Potentially Hazardeous Asteroids (PHA)
--  http://ssd.jpl.nasa.gov/sbdb_query.cgi?obj_group=pha;table_format=CSV;c_fields=AcBdBhBgBjBkBlBmApArAsAi;
--
--    ## for Near Earth Comets (NEC)
--  http://ssd.jpl.nasa.gov/sbdb_query.cgi?obj_group=neo;obj_kind=com;table_format=CSV;c_fields=AcBcBiBgBjBkBlBpApArAs;

-- 2)
--  The next thing to do is to press the "Generate table" button at the
--  end of the web page.
--  When prompted, choose "Save to disk" in your celestia root folder,
--  and GIVE IT THE APPROPRIATE NAME! ("asteroid.csv" or "comet.csv")

-- 3)
--  Run this script from Celestia. It should be in the Celestia/scripts
--  folder.

-- 4)
--  Once the script is done: quit Celestia and move the resulting files
--  from the Celestia root folder into your Celestia/extra folder

----------------------------------------------------------------------
-- NB:
-- The Period formula is different for hyperbolic vs. elliptic comets.
-- Comets with parabolic orbits will be turned into hyperbolic ones
-- by adding 0.00001 to their original eccentricity (=1).
-- See http://shatters.net/forum/viewtopic.php?t=8940 for details.

-- While some values are calculated, others, such as albedo and radius,
-- may sometimes only be estimated.
-- Therefore:
-- * Near Earth Asteroids without a known albedo or radius are given
--   a mean albedo of 0.11, which is then used to estimate the radius.
-- * Comets without a known radius are given an albedo of 0.03, and
--   their nuclei radius is estimated to an arbitrary (and underestimated)
--   average of 0.1km.

-- Nothing prevents you from modifying the script to run it with your
-- very own selection of data from the SBDB, but you should always keep
-- in mind that the files name need to be either "comet" or "asteroid",
-- that running the script on big files may produces a *HUGE* SSC file
-- (which is probably not yet suited for Celestia on our home computers...)
-- and, finally, that the average NEA albedo value of 0.11 will only
-- fit Near Earth objects, not farther objects.

-- The following pre-defined field sets are REQUIRED (and included in
-- the standard queries provided with this script):
--   ## For asteroids: "Keplerian elements" and "asteroid - physical"
--   ## For comets: "comet-style elements" and "comet - physical"

-- Enjoy!
--*********************************************************************


-------------------------------------------------
-- FUNCTIONS ------------------------------------
-------------------------------------------------


function celestia_keyboard_callback(input)
------------------------------------------
   key_pressed = input
   return(true)
end
---



function fromCSV (s, t)
-----------------------
-- Adapted from Roberto Ierusalimschy (LUA creator)
      s = s .. ','        -- ending comma
      local fieldstart = 1
      repeat
        -- next field is quoted? (start with `"'?)
        if string.find(s, '^"', fieldstart) then
          local a, c
          local i  = fieldstart
          repeat
            -- find closing quote
            a, i, c = string.find(s, '"("?)', i+1)
          until c ~= '"'    -- quote not followed by quote?
          if not i then error('unmatched "') end
          local f = string.sub(s, fieldstart+1, i-1)
          table.insert(t, (string.gsub(f, '""', '"')))
          fieldstart = string.find(s, ',', i) + 1
        else                -- unquoted; find next comma
          local nexti = string.find(s, ',', fieldstart)
          table.insert(t, string.sub(s, fieldstart, nexti-1))
          fieldstart = nexti + 1
        end
      until fieldstart > string.len(s)
 --     return t
end
---
 
 

function numF(leng, dec, n)
---------------------------
-- formats a number [n] for output, with a total
-- string length of [len] and [dec] digits float   
   
   return tostring(string.format("%"..leng.."."..dec.."f", tonumber(n)))
end
---



function remP(n)           
----------------
-- removes prepending punctuation caracters that
-- may interfere with math operations
   return string.gsub(n, '[<>=]*(.)', '%1')
end
---



function process_albedo_asteroid()
----------------------------------
   if (albedo == '') then
      albedo = 0.11
      albedo_message = "# a common average for Near Earth Asteroids"
   else
      albedo_message = "# from JPL SBDB"
   end
end
---


   
function process_albedo_comet()
-------------------------------
   if (albedo == '') then
      albedo = 0.03
      albedo_message = "# estimate of average albedo for active comets NUCLEI"
      -- From "Extreme albedo comets and impact hazard" (2004) W. M. Napier, J. T. Wickramasinghe and N. C. Wickramasinghe
   else
      albedo_message = "# from JPL SBDB"
   end
end
---
   

function process_diameter_asteroid()
------------------------------------
   if (diameter == '') then
         radius = 0.5 * 1329 *((math.pow(10, -(H)/5))/math.sqrt(albedo))
         radius_message = "# estimated: 0.5 * 1329 *(10^(-AbsMag[=H]/5))/sqrt(albedo))"
         -- http://www.physics.sfasu.edu/astro/asteroids/sizemagnitude.html
         -- http://cfa-www.harvard.edu/iau/lists/Sizes.html
   else
         radius = diameter/2
         radius_message = "# from JPL SBDB"
   end   
end
---


function process_diameter_comet()
---------------------------------
   if (diameter == '') then
         radius = 0.1
         radius_message = "# arbitrary constant (probably too small)"
         -- an arbitrary constant for comet nuclei radius
         -- (should be closer to 2km, but better undershoot than overshoot) 
   else
         radius = diameter/2
         radius_message = "# from JPL SBDB"
   end   
end
---



function process_e()
---------------------
   if (tonumber(e) < 1) then
        period = (math.pow(q/(1-e), 1.5))                     -- for elliptic orbits
        period_message = "# calculated: (q/(1-e))**1.5)"
   else
          if (e == 1) then e = 1 + 0.00001 end               -- turning parabolic orbits into hyperbolic ones (to avoid div by 0)
        period = (math.pow(q/(e-1), 1.5))                     -- for parabolic and hyperbolic orbits
        period_message = "# calculated: (q/(e-1))**1.5) (hyperbolic orbit)"
   end
end
---



function process_rot_per()
--------------------------
   if ((rot_per == '') or (tonumber(rot_per) == nil)) then
      rot_per = 1.2*math.exp(0.7+0.0012*(1+math.random(2300)))
--      random input into a formula approximating existing data
--      for 1970 asteroids.
      rot_per_message = "# random"
   else
      rot_per = remP(rot_per)
      rot_per_message = "# from JPL SBDB"
   end
end
---



function to_calendar_date(juliantimestamp)
------------------------------------------
   caltime = c:fromjulianday(tonumber(juliantimestamp))
   return string.format("%04d/%02d/%02d %02d:%02d:%02d", caltime.year, caltime.month, caltime.day, caltime.hour, caltime.minute, caltime.seconds)
end
---



function convert(input_name, input_extension)
---------------------------------------------
   os_date_stamp = os.date("*t")
   os_date_stamp = string.format("%04d_%02d_%02d", os_date_stamp.year, os_date_stamp.month, os_date_stamp.day)
   
   output_file = input_name.."_JPL_"..os_date_stamp..".ssc"
   input_file = input_name.."."..input_extension

   num = -1

 if (io.open(input_file,"r")) then
    io.close()
   output = io.open(output_file,"w")
   output:write("# Better check end of file for an incomplete object...")
   output:write("                                                         \r\n")
   -- a dummy line that will be erased later when the file is closed
   
   headers = {}

    if string.upper(input_name) == "COMET" then
      for line in io.lines(input_file) do
         COMET_print(line, input_name)
         c:print(header..message.." ...Processing "..input_file..": line "..num, 10, -1, 0)
         wait()
      end
    else
      for line in io.lines(input_file) do
         ASTEROID_print(line, input_name)
         c:print(header..message.." ...Processing "..input_file..": line "..num, 10, -1, 0)
         wait()
      end
    end
   c:print(header..message.."Cleaning file", 1, -1, 0)
   output:seek("set")
   output:write("# ",num," objects extracted from ",input_file," on ",string.gsub(os_date_stamp, "_", "/"),"         \r\n")
   output:flush()
   output:close()
   message = message.." "..num.."  objects written to  "..output_file.."\n"
 else
    message = message.." ****Oops! Couldn't find "..input_file..": "..output_file.." not created...\n"
 end
end
---



function ASTEROID_print(line, input_name)
-----------------------------------------


------------- Parameter gathering

if (num == -1) then
   fromCSV (line, headers)
   -- collect column headers
   num = num+1
else
   if line ~= '' then
      data = {}
      fromCSV (line, data)
      for k,v in pairs(data) do loadstring(headers[k]..' = "'..v..'"')() end
      -- collect each line from CSV file and transfer data into variables


------------- Data processing

      process_albedo_asteroid()
      process_diameter_asteroid()
      process_e()
      process_rot_per()
      

------------- Output      

      Designation = string.gsub(full_name, "^%s*(.-)%s*$", "%1")
      -- To remove spaces in front of the object's name
      Designation2 = string.gsub(Designation, "/", "")
      -- To avoid slashes in the object's name

         output:write("\n\n\n\n\n")
         output:write("# Asteroid ",Designation,"\n") -- Unaltered designation
   
         output:write("\"",Designation2,"\" \"Sol\"\n")
--         output:write("\"Ast",num,"\" \"Sol\"\n")
         output:write("   {\n")
         output:write("    Class \"asteroid\"\n")
         output:write("    Mesh \"asteroid.cms\"\n")
         output:write("    Texture \"asteroid.jpg\"\n")
         output:write("    EllipticalOrbit\n")
         output:write("       {\n")
         output:write("        Epoch             ",numF(16, 6, epoch),"\t\t# ",to_calendar_date(epoch),"\n")
         output:write("        Period            ",numF(25, 15, math.sqrt(math.pow(a,3))),"\t# calculated: sqrt(SemiMajorAxis**3)\n")
         output:write("        SemiMajorAxis     ",numF(25, 15, a)," \n")
         output:write("        Eccentricity      ",numF(25, 15, e)," \n")
         output:write("        Inclination       ",numF(25, 15, i)," \n")
         output:write("        AscendingNode     ",numF(25, 15, om)," \n")
         output:write("        ArgOfPericenter   ",numF(25, 15, w)," \n")
         output:write("        MeanAnomaly       ",numF(25, 15, ma)," \n")
         output:write("       }\n")
         output:write("    Radius                ",numF(12, 2, radius),"\t\t\t",radius_message,"\n")
         output:write("    Albedo                ",numF(12, 2, albedo),"\t\t\t",albedo_message,"\n")   
         output:write("    RotationPeriod        ",numF(12, 2, rot_per),"\t\t\t",rot_per_message,"\n")
         output:write("   }")

      num = num+1
      end
   end

end
---



function COMET_print(line, input_name)
--------------------------------------

--[[ From Selden Ball's page on how to transform orbital elements to Celestia's SSC formaat
   (http://www.lns.cornell.edu/~seb/celestia/transforming_ephemeris.html):
   
   "By convention, cometary elements don't include a mean anomaly value, so the Epoch you put into
   the SSC's EllipticalOrbit must be the time of the pericenter passage, when the mean anomaly is,
   by definition, 0.0."
]]

------------- Parameter gathering

if (num == -1) then
   fromCSV (line, headers)
   -- collect column headers
   num = num+1
else
   if line ~= '' then
      data = {}
      fromCSV (line, data)
      for k,v in pairs(data) do loadstring(headers[k]..' = "'..v..'"')() end
      -- collect each line from CSV file and transfer data into variables


------------- Data processing
   
      process_albedo_comet()
      process_diameter_comet()
      process_e()
      process_rot_per()
      M=0
      -- Mean anomaly at the time of perihelion

   
------------- Output
   
      Designation = string.gsub(full_name, "^%s*(.-)%s*$", "%1")
      -- To remove spaces in the object's name
      Designation2 = string.gsub(Designation, "/", " ")
      -- To replace slashes with spaces in the object's name

         output:write("\n\n\n\n\n")
         output:write("# Comet ",Designation,"\n") -- Unaltered designation
   
         output:write("\"",Designation2,"\" \"Sol\"\n")
--         output:write("\"Com",num,"\" \"Sol\"\n")
         output:write("   {\n")
         output:write("    Class \"comet\"\n")
         output:write("    Mesh \"halley.cmod\"\n")
         output:write("    Texture \"asteroid.jpg\"\n")
         output:write("    EllipticalOrbit\n")
         output:write("       {\n")
         output:write("        Epoch             ",numF(16, 6, tp),"\t\t# ",to_calendar_date(tp)," time of perihelion\n")
         output:write("        PericenterDistance",numF(25, 15, q),"\n")
         output:write("        Period            ",numF(25, 15, period),"\t",period_message,"\n")
         output:write("        Eccentricity      ",numF(25, 15, e)," \n")
         output:write("        Inclination       ",numF(25, 15, i)," \n")
         output:write("        AscendingNode     ",numF(25, 15, om)," \n")
         output:write("        ArgOfPericenter   ",numF(25, 15, w)," \n")
         output:write("        MeanAnomaly       ",numF(9, 0, M)," \n")
         output:write("       }\n")
         output:write("    Radius                ",numF(12, 2, radius),"\t\t\t",radius_message,"\n")
         output:write("    Albedo                ",numF(12, 2, albedo),"\t\t\t",albedo_message,"\n")   
         output:write("    RotationPeriod        ",numF(12, 2, rot_per),"\t\t\t",rot_per_message,"\n")
         output:write("   }")

      num = num+1
   end
end

end
---



-------------------------------------------------
-- SCRIPT START----------------------------------
-------------------------------------------------
script_version = "JPL SBDB .csv to Celestia .ssc converter (beta1)"
header = " "..script_version.."\n __________________________________________\n"
c = celestia
peri = 0

c:requestkeyboard(true)
while(key_pressed ~= " ") do
   c:print(header.." This script will look for .CSV files in your\n Celestia root folder and convert them to\n Celestia .SSC files.\n\n You will need to create and download\n at least one .csv file (comet.csv or asteroid.csv)\n using the JPL SBDB query page:\n http://ssd.jpl.nasa.gov/sbdb_query.cgi\n\n If you already downloaded one or both files,\n press the Spacebar to run the conversion.\n Otherwise, press the Escape key to exit now...", 1, -1, 0)
   wait()
end
c:requestkeyboard(false)
c:print(header.." Celestia will now request system access\n in order to handle the output file(s).", 5, -1, 0)
wait(4)
c:requestsystemaccess()
wait()
num = 0
message = ""
input_extension = ""




---------------------------------
-- CONVERSION -------------------
---------------------------------
headers = {}        -- table to collect headers
data = {}           -- table to collect fields
-- the data file should resemble the following:
--line1 = 'full_name,name,prefix,pha,full_name,e,q,i,om,w,ad,tp_cal,per_y,class,data_arc,condition_code,n_obs_used,two_body,A1,A2,A3,DT'
--line>1 = '"     C/-146 P1",,C,,"     C/-146 P1",1,0.43,71,330,261,,-1460628.0000000,,PAR,5,,4,T,,,,'


convert('comet', 'csv')
convert('asteroid', 'csv')

c:requestkeyboard(true)
key_pressed = nil
while(key_pressed == nil) do
   c:print(header..message.."\n Script complete, press any key to return to Celestia", 1, -1, 0)
   wait()
end
c:requestkeyboard(false)


--[[
----------------------------
Typical .ssc asteroid item :
----------------------------
"Ceres" "Sol"
   {
    Class "asteroid"
   Texture "asteroid.jpg"
   Radius 487.5        # 975x909km
   Oblateness 0.068    # Thomas et al, Nature, 8 Sep 2005
   EllipticalOrbit
      {
       Epoch       2452600.5     # 2002 Nov 22 00:00UT
      Period            4.60014
      SemiMajorAxis     2.7660
      Eccentricity      0.0793
      Inclination      10.584
      AscendingNode    80.483
      ArgOfPericenter  74.043
      MeanAnomaly     232.067
      }
   RotationPeriod         9.075
   RotationEpoch    2449249.91125 # 1993 Sep 19 09:52:12UT
   Obliquity             11
   EquatorAscendingNode  29
   RotationOffset       339.85
   Albedo 0.113
   }



Typical .ssc comet item :
----------------------------
"C1995 O1 (Hale-Bopp)" "Sol"
{
   Class "comet"
   Mesh "asteroid.cms"
   Texture "asteroid.jpg"
   Radius 50.0      # guess
        Albedo 0.1
   EllipticalOrbit
   {
      Period             2550.259099
      PericenterDistance    0.931444
      Eccentricity          0.995010
      Inclination          89.36570
      AscendingNode       282.50470
      ArgOfPericenter     130.85030
    MeanAnomaly           0.0
      Epoch           2450537.8292
   }
   InfoURL "http://ssd.jpl.nasa.gov/cgi-bin/eph"
}

]]