src/igrf2c.py

   1 #!/usr/bin/python3
   2 # Martin Budaj 2014
   3
   4 import sys
   5
   6 if len(sys.argv) == 1:
   7   print('usage: %s <coefficients.txt>' % sys.argv[0])
   8   sys.exit(1)
   9
  10 def get_values(d):
  11   return d[3:len(d)-1]
  12
  13 def iround(d):
  14   if '.' in d:
  15     return d.rstrip('0').rstrip('.')
  16   else:
  17     return d
  18
  19 data = dict()
  20 max_n = 0
  21
  22 with open(sys.argv[1]) as fin:
  23   for l in fin:
  24     if l.startswith(('#', 'c/s')): continue
  25     l = l.rstrip()
  26     ldata = l.split()
  27     if ldata[0] == 'g/h':
  28       assert(ldata[1] == 'n' and ldata[2] == 'm')
  29       years = get_values(ldata)
  30     elif ldata[0] in ('g', 'h'):
  31       n = int(ldata[1])
  32       m = int(ldata[2])
  33       data[ldata[0],n,m] = list(get_values(ldata))
  34       data[ldata[0]+'_delta',n,m] = ldata[len(ldata)-1],
  35       max_n = max(max_n, n, m)
  36     else:
  37       raise ValueError('data line in unknown format!')
  38
  39 with open('../thgeomagdata.h','w') as fout:
  40   fout.write('''// generated by geomag/igrf2c.py
  41
  42 #ifndef thgeomagdata_h
  43 #define thgeomagdata_h
  44
  45 ''')
  46   for gh in ('g', 'h', 'g_delta', 'h_delta'):
  47     if 'delta' in gh:
  48       isdelta = True
  49       suffix = 'D'
  50       yy = (0,)
  51     else:
  52       isdelta = False
  53       suffix = '[%d]' % len(years)
  54       yy = range(0,len(years))
  55     fout.write('static const double thgeomag_%sNM%s[%d][%d] = {\n' %
  56       (gh[0].upper(),suffix,max_n+1,max_n+1))
  57     for y in yy:
  58       if not isdelta: fout.write('{\n')
  59       for n in range(0,max_n+1):
  60         fout.write('  {')
  61         for m in range(0,max_n+1):
  62           fout.write('%s, ' % (iround(data[gh,n,m][y]) if (gh,n,m) in data else '0'))
  63         fout.write('},\n')
  64       if not isdelta: fout.write('},\n')
  65     fout.write('};\n')
  66
  67   fout.write(
  68 '''#define thgeomag_maxmindex %d
  69 #define thgeomag_step %d
  70 #define thgeomag_minyear %d
  71 #define thgeomag_maxdeg %d
  72
  73 #endif
  74
  75 ''' % (len(years)-1, float(years[1])-float(years[0]), float(years[0]), max_n))
  76
  77   print('OK: ../thgeomagdata.h created')