#!/usr/local/bin/python

import os, sys
import math

import Numeric as Nu

import nast.gomodelRNAmaj as gomodel
import nast.cubicSplineSetup as cubicSplineSetup

def cubicSplineCalc(xa, ya, y2a, x):
  n = len(xa)
  klo=0
  khi=n-1
  while (khi-klo > 1):
    k = (khi+klo)/2
    if (xa[k] > x): khi=k
    else: klo=k

  h=float(xa[khi]-xa[klo])
  if (h == 0.0): return 1
  hScale=1/h

  a = hScale*(xa[khi]-x)
  dadx =-hScale

  b = hScale*(x-xa[klo])
  dbdx = hScale

  y  = a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0
  yp = dadx*ya[klo]+dbdx*ya[khi]+((3*a*a*dadx-dadx)*y2a[klo]+(3*b*b*dbdx-dbdx)*y2a[khi])*(h*h)/6.0

  return (y, yp)

def usageError():
    print 'usage: %s splineInFilename splineOutFilename \\' \
         % os.path.basename(sys.argv[0])
    print '                                  [splineFirstDerivativeOutFile]'
    sys.exit(1)


if __name__=='__main__':
    if len(sys.argv)<3: usageError()
    splineFilename = sys.argv[1]
    plotFilename = sys.argv[2]
    fOut = open(plotFilename, 'w')

    try:
        plotFilename2 = sys.argv[3]
        fOut2 = open(plotFilename2, 'w')
    except IndexError:
        plotFilename2 = ''
  
    splineX, splineY, splineIndex = gomodel.readSplineFiles([splineFilename])[0]
    splineY2= cubicSplineSetup.makeY2(splineX, splineY)

    xypSpline = []
    numPoints = 100
    for i in range(numPoints):
        x = i*(math.pi/(numPoints-1))
        y, yp = cubicSplineCalc(splineX, splineY, splineY2, math.cos(math.pi-x))
        fOut.write('%.3f  %.3f\n' % (x, y))
        if plotFilename2:
            fOut2.write('%.3f  %.3f\n' % (x, yp/20))

    fOut.close()