#!/usr/bin/python
# -*- coding: UTF-8 no BOM -*-

import os,string,sys,re
from optparse import OptionParser
import damask

scriptID   = string.replace('$Id$','\n','\\n')
scriptName = scriptID.split()[1]

sampleSym = {  'Orthotropic'  : (90,90,90),
               'Triclinic'    : (360,180,360)
            }

# --------------------------------------------------------------------
#                                MAIN
# --------------------------------------------------------------------

parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Transform the OIM texture data into linear ODF data, the input file can be generated by software "TSL OIM Analysis",
and the output file can be used by the script: "hybridIA_linODFsampling.py"
""", version=string.replace(scriptID,'\n','\\n')
)

parser.add_option('-f', '--file', dest='file', type='string', metavar = 'string', \
                     help='file name')
parser.add_option('-s', '--symmetry', dest='symmetry', choices=sampleSym.keys(),\
                     help='Sample symmetry, the default is [%default]')
parser.add_option('--step', dest='step', type='int', nargs=3,
                     help='''Angle steps of: phi1, Phi, phi2. For example, if the 
                           range of phi1 is 360 degree, and the desired interval 
                           is 5 degree, then the angle step of phi1 is 72, but 
                           we will discard the 72th step, because 360 is 
                           identical with 0, the default depends on the sample
                           symmetry. ''', \
                           metavar='int int int')
                                     
parser.set_defaults(symmetry = 'Triclinic')
parser.set_defaults(step     = None)
options = parser.parse_args()[0]

if options.step is None:
  options.step= [sampleSym[options.symmetry][i]/5 for i in xrange(3)]

if not os.path.exists(options.file):
  parser.error('texture file does not exist'); sys.exit()

inName = options.file
outName = os.path.splitext(inName)[0]+'.linearODF'
nPhi1,nPHI,nPhi2 = options.step
dPhi1,dPHI,dPhi2 = [sampleSym[options.symmetry][i]/options.step[i] for i in xrange(3)]

N = (nPhi1-1)*(nPHI-1)*(nPhi2-1)

try:
  inFile = open(inName,'r')
  content = inFile.readlines()
except:
  print 'unable to read:',inName
  sys.exit(1)
try:
  outFile = open(outName,'w')
except:
  print 'unable to write:',outName
  sys.exit(1)

ODF = [[[[None] for k in range(nPhi2)] for j in range(nPHI)] for i in range(nPhi1)]
linear = [None]*N
line = 0

while (content[line].startswith('#')):  # skip comments at start of file
  line += 1

for iPhi1 in range(nPhi1):
  for iPHI in range(nPHI):
    for iPhi2 in range(nPhi2):
      words = content[line].split()
      ODF[iPhi1][iPHI][iPhi2] = float(words[3]) # extract intensity (in column 4)
      line += 1

for iPhi1 in range(nPhi1-1):
  for iPHI in range(nPHI-1):
    for iPhi2 in range(nPhi2-1):
      linear[iPhi1*(nPHI-1)*(nPhi2-1)+iPHI*(nPhi2-1)+iPhi2] = (\
        ODF[iPhi1  ][iPHI  ][iPhi2  ] + \
        ODF[iPhi1  ][iPHI  ][iPhi2+1] + \
        ODF[iPhi1  ][iPHI+1][iPhi2  ] + \
        ODF[iPhi1  ][iPHI+1][iPhi2+1] + \
        ODF[iPhi1+1][iPHI  ][iPhi2  ] + \
        ODF[iPhi1+1][iPHI  ][iPhi2+1] + \
        ODF[iPhi1+1][iPHI+1][iPhi2  ] + \
        ODF[iPhi1+1][iPHI+1][iPhi2+1] \
      ) / 8.0


inFile.close()
outFile.write('%-6i%-6i%-6i  #Ranges of phi1, Phi, phi2\n'%sampleSym[options.symmetry])
outFile.write('%-6.2f%-6.2f%-6.2f  #Deltas of phi1, Phi, phi2\n'%(dPhi1,dPHI,dPhi2))
outFile.write('%-6i%-6i%-6i  #Angular steps needed to be converted\n'%(nPhi1-1,nPHI-1,nPhi2-1))
outFile.write('cell-centered data\n')
outFile.write('\n')

for i in range(N):
  outFile.write('%g\n'%(linear[i]))
  
outFile.close()