#!/usr/bin/env python
# -*- encoding:utf-8 -*-

import math

class MyLine:
  def __init__( self, gr = 0.0, inter = 0.0 ):
    # y = gradient * x + intercept
    self.gradient = gr
    self.intercept = inter

  def getValue( self, x ):
    return self.gradient * x + self.intercept

  def crosspoint( self, myl ):
    if not isinstance( myl, MyLine ):
      return
    return self.crosspoint_extended( myl.gradient, myl.intercept )

  def crosspoint_extended( self, gr, inter ):
    del_grad = self.gradient - gr
    del_inter = inter - self.intercept
    if del_grad == 0.0:
      return
    crossx = del_inter / del_grad
    crossy = self.getValue( crossx )
    return [ crossx, crossy ]

  # generate MyLine which passes (px,py) and vertical to this Line
  def genLineFromPoint( self, px, py ):
    ret = MyLine()
    ret.gradient = -1.0 / self.gradient
    ret.intercept = py - ret.gradient * px
    return ret


def unroll_angle( angle ):
  angle -= 360 * round( angle / 360 )
  return angle

if __name__ == "__main__":
  # alpha helix
  equation = MyLine( -1, -105 )
  chara_len = 45
  phi_range = [ -90, -35 ]
  ## 3-10 helix
  #equation = MyLine( -1, -75 )
  #chara_len = 50
  #phi_range = [ -90, 15 ]
  ## pi helix
  #equation = MyLine( -1, -125 )
  #chara_len = 40
  #phi_range = [ -85, -40 ]
  ## extended beta
  #equation = MyLine( -1, 0 )
  #chara_len = 70
  #phi_range = [ -140, -100 ]
  ##### common #####
  psi_range = [ equation.getValue( phi_range[0] ), equation.getValue( phi_range[1] ) ]
  alpha = 1.0 / ( chara_len * chara_len )
  phis = []
  psis = []
  for i in range( -18, 18 ):
    phis.append( 10 * i )
    psis.append( 10 * i )

  for phi in phis:
    for psi in psis:
      myl = equation.genLineFromPoint( phi, psi )
      cp = equation.crosspoint( myl )
      # distance from the crosspoint
      delphi = unroll_angle( phi - cp[0] )
      delpsi = unroll_angle( psi - cp[1] )
      distsq = delphi * delphi + delpsi * delpsi
      if unroll_angle( cp[0] - phi_range[0] ) < 0:
        delphi2 = unroll_angle( cp[0] - phi_range[0] )
        delpsi2 = unroll_angle( cp[1] - psi_range[0] )
        distsq += delphi2 * delphi2 + delpsi2 * delpsi2
      elif unroll_angle( cp[0] - phi_range[1] ) > 0:
        delphi2 = unroll_angle( cp[0] - phi_range[1] )
        delpsi2 = unroll_angle( cp[1] - psi_range[1] )
        distsq += delphi2 * delphi2 + delpsi2 * delpsi2
      print phi, psi, math.exp( - alpha * distsq )
    print
