# 
# calculate duple orientations by RDC periodicity

import os
import numpy
from math import radians, degrees,sin, cos, pi, sqrt 
from lib.levmar import dbcdif
ALF_PM = pi*2.0/11.0
pi2 = pi*2 
from random import random 
from libio import MIO
from copy import copy 
from librdc import DuplexRDC

def rdcwave(parameters, dummy, ndata):
    """ Da, R, Theta, Phi, Rho are unknow """
    Da,  R, the, phi, rho_0, delta = parameters
    pm = the, phi, rho_0, delta, Da, R, ndata
    return pm2rdcs_wave(pm)

 
def pm2rdcs_wave(parameters):
    the, phi, rho_0, delta, Da,R, ndata = parameters
    sin_del = sin(delta)
    cos_del = cos(delta)
    sin_2del = 2*sin_del*cos_del
    Ki = (3.0 /16.) * (4.0 + 6.0*R*cos(2.*phi) + R*cos(2.*(phi-the)) - 4. *cos(2.*the) + R*cos(2.*(phi+the)) )*sin_del*sin_del 
    Kii = - Ki
    Kiii = -3.0 * R * cos(the)*sin(2. * phi)* sin_del * sin_del
    Kiv =  0.75  * (-2. + R*cos(2.* phi) )*sin(2.*the)* sin_2del
    Kv = - 3.0 * R* sin(phi)*sin(the)*cos(phi)* sin_2del
    Kvi = (1.0/32.) * ( 4. + 6.* R* cos(2.*phi) - 3. *R*cos(2.*(phi-the)) + 12.*cos(2.*the) - 3.*R*cos(2.*(phi+the)) )*(1. + 3.*cos(2.*delta))
    rdcs = []
    cst =  Kvi + 0.5*(Ki + Kii) 
    for n in ndata:#[1:]
        alf = ALF_PM * (n-1)
        rho = rho_0 + alf 
        rdc  = Da*( 0.5 * Kiii * sin(2.0*rho) + 0.5* (Ki - Kii)*cos(2.*rho) + Kv * sin(rho) + Kiv*cos(rho) + cst )
        rdcs.append(rdc)
    return rdcs 

def SIMRDCWave(parameters, dummy, ndata):
    num_dps = int((len(parameters) -2 )/ 4 ) #
    Da = parameters[-2]
    R = parameters[-1]
    cnt = 0
    nd_id = 0
    rdcs = []
    #print "# all_ndata: ", ndata
    for i in range(num_dps):
        the = parameters[cnt]
        phi = parameters[cnt+1]
        rho = parameters[cnt+2]
        delta = parameters[cnt+3]
        cnt += 4
        num_nd = ndata[nd_id]
        #print "# num_nd: ", num_nd
        nd_id += 1
        nd_end = nd_id + num_nd
        nd = ndata[nd_id:nd_end]
        nd_id = nd_end
        pm = the, phi, rho, delta, Da, R, nd 
        rdcs += pm2rdcs_wave(pm)
        #print "# pm: ", pm 
    #print rdcs 
    return rdcs 


class RDCFit(MIO):
    """ fit rdc by rdc periodicity """
    fit_full_keys  = ["Da", "R", "phi", "the", "rho","delta"]
    def bc2IFitVar(self, lb, ub):
        """ return inital random values of the variables in their ranges for fitting """
        num_var = len(lb)
        x0 = [0]*num_var
        for i in range(num_var):
            assert ub[i] > lb[i]
            x0[i] = lb[i] + (ub[i]-lb[i]) * random()
        return x0 

    
    def getLU(self, fitpm, keys=None, vb=0):
        """ self.parameters = [ Da, R, the, phi, rho_0] """
        lb = []
        ub = []
        if keys is None:
            keys = ['Da', 'R','the','phi', 'rho', 'delta']
        for _key in keys: 
            if _key in fitpm:
                _val = fitpm[_key]
                lb.append(_val[0])
                ub.append(_val[1])
                if vb: print "# %5s: min= %7.3f  max= %7.3f " %  (_key, degrees(_val[0]), degrees(_val[1]))
        return lb, ub 
    
    
    def sortByRmsd(self, rmsds, all_results):
        #print "rmsds: ", rmsds 
        rmsds = numpy.array(rmsds)
        sids = numpy.argsort(rmsds)        
        num = min(50, len(rmsds))
        #print "# num= %3d " % num 
        for i in range(num):
            sid = sids[i]
            _result = all_results[sid][0]
            #print "# sid= %3d " % sid 
            rmsd = rmsds[sid]
            self.prnFitResults(_result, rmsd, rho_chk=0)
        return all_results[sids[0]] #$[0]                  
        
            
    
    def chkVeglia(self, parameters, vb=0):
        rdcs = rdcwave(parameters,None,self.ndata)
        assert len(self.ydata) == len(rdcs)
        num_data = len(self.ydata)
        sum_sq = 0.
        lns = []
        rdc_difs = []
        for i in range(num_data):
            diff =  self.ydata[i]-rdcs[i]
            rdc_difs.append(diff)
            sum_sq += diff * diff
            lns += [ "# bpid= %3d resid= %3d  exp_rdc= %8.3f  cal_rdc= %8.3f  df= %8.3f " % (self.bp_ids[i], self.dp_rdc_resids[i], self.ydata[i], rdcs[i], self.ydata[i]-rdcs[i])]
        rmsd = sqrt(sum_sq/float(num_data))
        lns += [ "# rmsd= %8.3f,  Da= %8.5f  R= %8.4f   the= %6.1f  phi= %6.1f rho= %6.1f delta= %6.1f  " % \
                 (rmsd, parameters[0], parameters[1], degrees(parameters[2]), degrees(parameters[3]), degrees(parameters[4]), degrees(parameters[5]))]
        if vb: self.prnLns(lns)
        return rmsd, sum_sq, rdc_difs, lns  
    
    
    def fitRDC(self):
        """
        parameters = [ Da, R, the, phi, rho_0]
        """
        lb, ub = self.getLU(self.fitpm)
        all_results = []
        rmsds = []
        for i in range(self.try_num):
            x0 = self.bc2IFitVar(lb,ub)
            results = dbcdif(rdcwave, x0, self.ydata, lb,ub, self.max_iter, data=self.ndata)
            rmsd, sum_sq, rdc_difs, lns = self.chkVeglia(results[0], vb=0)
            #self.prnFitResults(results, rmsd, rho_chk=0)
            all_results.append((results,rdc_difs, rmsd, lns))
            rmsds.append(rmsd)
        return self.sortByRmsd(rmsds, all_results)
    
    def prnFitPM(self):
        """ print out RDC periodicity fitting parameters """
        lns = ["# RDC periodicity fitting parameters "]
        fit_pm = self.fitpm
        pm_keys = fit_pm.keys()
        for _key in self.fit_full_keys:
            if _key in pm_keys:
                if _key in ['Da', "R"]:
                    lns += ["# FitPM:  %-5s= (%7.2f -- %7.2f) " % (_key, fit_pm[_key][0], fit_pm[_key][1]) ]
                else:
                    lns += ["# FitPM:  %-5s= (%7.1f -- %7.1f) " % (_key, degrees(fit_pm[_key][0]), degrees(fit_pm[_key][1])) ]
        lns += ["# fitting try number= %3d  maximum_iteration= %3d " % (self.try_num, self.max_iter) ]
        self.add2log(lns)
        
        
    def backChkAxis(self, parameters, helix=2):
        #parameters =   Da,  R, the, phi, rho_0, delta = parameters
        print "\n# back calculating RDCs "
        #print "# pm: ", parameters
        #rmsd, sum_sq,  lns = self.chkVeglia(parameters, vb=1)
        self.prnLns(parameters[-1])
        #self.estimateRDC(parameters, self.ndata)

    def estimateRDC(self, pm, ndata):
        bpids = ndata[2:]
        
        full_bpids = range(bpids[0], bpids[-1]+1)
        #full_bpids = range(1, bpids[-1]+1)
        est_bpids = []
        for bpid in full_bpids:
            if bpid not in bpids:
                est_bpids.append(bpid)
        num_est_bpids = len(est_bpids)
        if num_est_bpids > 0:
            new_ndata = ndata[:2] + est_bpids
            rdcs = veglia(pm,None, new_ndata)
            print "\n# estimated RDCs by Duplex RDC periodicity "
            print "# Da= %8.3f  R= %8.3f the= %6.1f phi= %6.1f rho= %6.1f delta= %6.1f " % \
                  (pm[0], pm[1],  degrees(pm[2]), degrees(pm[3]), degrees(pm[4]), degrees(pm[5]))
            for i in range(num_est_bpids):
                print "#  bpid= %3d  cal_rdc= %8.3f "  % (est_bpids[i], rdcs[i])
            
    def prnFitResults(self, results, rmsd, rho_chk=1):
        #print results 
        Da,R,the,phi,rho_0, delta  = results[0]
        if rho_chk:
            if degrees(rho_0) < 180:
                print "# rmsd= %8.3f  Da= %8.3f  R= %8.3f "  % (rmsd, Da, R), 
                print "# the= %6.1f  phi= %6.1f  rho_= %6.1f  delta= %6.1f" % (degrees(the), degrees(phi), degrees(rho_0), degrees(delta))
        else:
            print "# rmsd= %8.3f  Da= %8.3f  R= %8.3f "  % (rmsd, Da, R), 
            print "# the= %6.1f  phi= %6.1f  rho_= %6.1f  delta= %6.1f" % (degrees(the), degrees(phi), degrees(rho_0), degrees(delta))            
        
    
    def setFitDa_Range(self, err_range=(-30, 30)):
        self.Da_err = err_range 
            
    def setFitR_Range(self, err_range=(0, 0.667)):
        self.R_err = err_range 
    
    def setFitPhi_Range(self, err_range=(0, pi)):
        self.Phi_err = err_range 
            
    def setFitThe_Range(self, err_range=(0, pi/2.0)):
        self.The_err= err_range
            
    def setFitRho_Range(self, err_range=(0, pi2)):
        self.Rho_err= err_range
        
    def setFitDelta_Range(self, err_range=(radians(97), radians(108))):
        self.Delta_err= err_range    
        #keys=["Da", "R", "phi", "the", "rho","delta"]
        
    def setFitErr_range(self):
        """ set default error range for RDC periodicity fitting """
        self.setFitDa_Range()
        self.setFitR_Range()
        self.setFitPhi_Range()
        self.setFitThe_Range()
        self.setFitRho_Range()
        self.setFitDelta_Range()
    
    
    def setFitPM_DaR(self, Da=None, R=None, try_num=20, max_iter=200):
        if Da is None:
            print "# must in Da and R values "
            raise 
        self.Da  = Da 
        self.R = R 
        keys = ['the', 'phi','rho','delta']
        self.setFitPM(keys=keys,try_num=try_num, max_iter=max_iter)
            
        
    def setFitPM(self, keys= None, try_num=20, max_iter=200):    
        fitpm = {}
        if keys is None: keys = self.fit_full_keys
        for _key in keys:
            if _key == "Da":
                fitpm[_key] = self.Da_err
            elif _key == "R":
                fitpm[_key] = self.R_err
            elif _key == "phi":
                fitpm[_key] = self.Phi_err
            elif _key == 'rho':
                fitpm[_key] = self.Rho_err
            elif _key == 'delta':
                fitpm[_key] = self.Delta_err
            elif _key == 'the':
                fitpm[_key] = self.The_err
        self.fitpm = fitpm
        self.setFitYdata(keys=keys)
        self.setTryIterNum(try_num=try_num, max_iter=max_iter)
    
    def setFitYdata(self, keys=None):
        if keys is None: keys = self.fit_full_keys
        self.ndata = copy(self.bp_ids)
        self.ydata = copy(self.bp_rdcs)
        if "Da" not in keys: self.ndata = [self.Da, self.R] + self.ndata 
        
        
    def setTryIterNum(self, try_num=20, max_iter=200):
        self.try_num=try_num
        self.max_iter = max_iter
         
    
class RDCSIMFIT(RDCFit, DuplexRDC):
    """ simultaneously fit all duplexes' RDCs """
    max_iter = 200
    def setSimBPIds(self):
        """ set all duplex'es basepair IDs for rdc periodicty fitting """
        self.all_rdcs = []
        self.all_rdcs_resids = []
        self.all_bp_ids = []
        self.ndata = []
        self.all_dup_presids = {} # all duplexs' phase resids 
        
        self.dups_rdcs = {}
        self.dups_bpids = {}
        self.dups_rdc_resids = {}
        
        for dup_id in  self.duplex_ids:
            self.setBpIds(dup_id, vb=1)
            #print self.bp_ids
            #print self.bp_rdcs
            #print self.dp_rdc_resids
            self.all_dup_presids[dup_id] = self.dup_phase_resid
            self.all_rdcs +=  self.bp_rdcs
            self.dups_rdcs[dup_id] = self.bp_rdcs
            
            self.all_rdcs_resids += self.dp_rdc_resids
            self.dups_rdc_resids[dup_id] = self.dp_rdc_resids 
            
            self.dups_bpids[dup_id] = self.bp_ids
            self.ndata += [len(self.bp_ids)] + self.bp_ids
            self.all_bp_ids += self.bp_ids
        #print self.all_bp_ids
          
    
    def setSimFitPM(self, TEST=1):
        fit_pm = []
        pi2 = pi*2.0 
        for dup_id in self.duplex_ids:
            if dup_id == '1' and TEST: 
                #print "####DD1"
                the1, the2 = 20, 28
                phi1, phi2 = 85, 100
                fit_pm.append((radians(the1), radians(the2))) #print "dd1"
                fit_pm.append((radians(phi1), radians(phi2)))
                lns = [ "# duplex %s: the: (%6.1f--%6.1f) phi: (%6.1f--%6.1f) " % (dup_id,the1, the2, phi1, phi2)]
                self.prnLines(lns)
                self.add2log(lns)
            elif dup_id == '2' and TEST:
                print "# ddd2::"
                the1, the2 = 5,15
                phi1, phi2 = 135, 145
                fit_pm.append((radians(the1), radians(the2))) #print "dd1"
                fit_pm.append((radians(phi1), radians(phi2)))
                lns = [ "# duplex %s: the: (%6.1f--%6.1f) phi: (%6.1f--%6.1f) " % (dup_id,the1, the2, phi1, phi2)]
                self.prnLines(lns)
                self.add2log(lns)
            elif dup_id == '3' and TEST: 
                print "# ddd3::"
                the1, the2 = 0, 10
                phi1, phi2 = 30, 40
                fit_pm.append((radians(the1), radians(the2))) #print "dd1"
                fit_pm.append((radians(phi1), radians(phi2)))
                lns = [ "# duplex %s: the: (%6.1f--%6.1f) phi: (%6.1f--%6.1f) " % (dup_id,the1, the2, phi1, phi2)]
                self.prnLines(lns)
                self.add2log(lns)
            else:    
                fit_pm.append((0, pi/2.0)) # the 
                fit_pm.append((0, pi))  # phi
            fit_pm.append((0, pi2)) 
            fit_pm.append((radians(97), radians(108))) 
            #fit_pm.append((radians(101), radians(105))) 
        fit_pm.append((-30.0, 30))      # Da
        fit_pm.append((0.001, 0.667))   # R 
        self.fit_pm = fit_pm 
    
    
    def simPM2LU(self):
        """ self.parameters = [ Da, R, the, phi, rho_0] """
        lb = []
        ub = []
        for _ll, _uu in self.fit_pm:
            lb.append(_ll)
            ub.append(_uu)
        return lb, ub 
    
    def simTest(self):
        lb, ub = self.simPM2LU() #getLU(fitpm, keys=pm_keys, vb=1)
        #all_results = []
        #rmsds = []
        #for i in range(self.num_try):
        x0 = self.bc2IFitVar(lb,ub)
        SIMRDCWave(x0, None, self.all_bp_ids)
        
        
    def simFitRDC(self, vb=0):
        """
        """
        
        lb, ub = self.simPM2LU() 
        all_results = []
        rmsds = []
        lns = ["# Levmar fit iteration: %4d   max_num_try= %4d " % (self.max_iter, self.num_try) ]
        self.add2log(lns, vb=1)
        print "# runing simultaneous RDC fitting ... take a long time "
        for i in range(self.num_try):
            print "# Fitting iteration:  %4d/%d " % (i+1, self.num_try)
            x0 = self.bc2IFitVar(lb,ub)
            results = dbcdif(SIMRDCWave, x0, self.all_rdcs, lb,ub, self.max_iter, data= self.ndata)
            rmsd, sum_sq,  lns = self.chkSimFitRes(results[0], vb=0)
            self.setSimDupAxle(results[0])
            if not self.ANGFILT or self.simFiltDupAng():
                angs_dup = self.getSimDupAng()
                all_results.append((results, angs_dup, lns))
                rmsds.append(rmsd)
                #if rmsd < 2.0:
                #    self.add2log(lns)
        self.sim_fit_res =  self.sortSimFitByRmsd(rmsds, all_results, vb=vb)
    
    
    def saveSimFitRes(self, fn_pre=None):
        """ self.sim_res = [(rmsd, pm, angs_dup), ...] """
        if fn_pre is None:
            fn_pre = self.fn_pre 
        fn = "%s_simfit.bin" % fn_pre
        self.dumpObj(self.sim_fit_res, fn)
    
    def prnSimFitRes(self, sim_fit_res, prn=1):
        rmsd, pm, angs_dup, rdc_lns  = sim_fit_res
        Da,R,tprs = self.simPM2DPR_DaR(pm)
        lns =  ['## Best simultaneous fitting result ']
        lns += ['## RDC rmsd= %8.3f   Da= %8.3f  R= %6.2f ' % (rmsd,Da,R) ]
        lns += ['# duplex Orientations from Orient module']
        lns += ['## dup_id  the    phi   rho    delta phase_resid']
        for i in range(len(self.duplex_ids)):
            dup_id = self.duplex_ids[i]
            phase_resid = self.all_dup_presids[dup_id]
            tpr_deg = [ degrees(_tpr) for _tpr in tprs[i]]
            lns += ["%7s %6.0f %6.0f %6.f %6.0f  %4d "  % (dup_id, tpr_deg[0], tpr_deg[1], tpr_deg[2], tpr_deg[3], phase_resid) ]
        if prn: self.prnLines(lns)
        self.add2log(rdc_lns)
        return lns 
        
    def saveBestFit(self, vb=0):
        """ self.sim_fit_res = [(rmsd, pm, angs_dup), ...] """
        lns = self.prnSimFitRes(self.sim_fit_res[0], prn=0)
        self.add2log(lns)
        if vb: self.prnLines(lns)
        self.saveLines(lns, "%s_best_sim_fit.txt" % self.fn_pre, vb=vb)
           
    
    def sortSimFitByRmsd(self, rmsds, all_results, **rest):
        if len(rmsds)  == 0: self.err_exit("# could not find accepted fitting results: increase rmsd cutoff or angle range ")
        rmsds = numpy.array(rmsds)
        sids = numpy.argsort(rmsds)        
        if not hasattr(self, 'num_saved_fit'):
            self.num_saved_fit = 100
        num = min(self.num_saved_fit, len(rmsds))
        best_res = []
        best_rmsd_lns = []
        for i in range(num):
            sid = sids[i]
            result, angs_dup, lns = all_results[sid]
            pm = result[0]
            rmsd = rmsds[sid]
            best_rmsd_lns.append(self.prnSimFitResult(pm, rmsd, angs_dup, **rest))
            self.add2log(lns)
            best_res.append((rmsd, pm, angs_dup, lns))
        self.add2log(best_rmsd_lns)
        return best_res
        
    
    def prnSimFitResult(self, pm, rmsd, ang_dups, **rest):
        ln = " %7.3f : " % rmsd
        fmt = "%6.0f " * len(ang_dups)
        ln +=  fmt % ang_dups
        fmt = ": " +  "%5.0f " * (len(pm) - 2)
        pm_deg = [ degrees(_tpr) for _tpr in pm[:-2]]
        ln  += fmt % tuple(pm_deg)
        ln +=  ": %6.1f %5.2f" % (pm[-2], pm[-1])
        vb = rest.get('vb', 0)
        if vb: print ln 
        self.add2log([ln])
        return ln 
    
    
    def thePhi2Axis(self, the, phi):
        """ the, phi in units radians """
        #the = radians(the)
        #phi = radians(phi)
        #print "# the:", the, "phi: ", phi 
        x = cos(phi)*sin(the)
        y = sin(phi) * sin(the)
        z = cos(the)
        return numpy.array([x,y,z])
    
            
    def chkSimFitRes(self, pms, vb=0):
        rdcs = SIMRDCWave(pms, None,self.ndata)
        exp_rdcs = self.all_rdcs
        num_rdc = len(exp_rdcs)
        #if 1: print  "# num_rdcs in Simfit: %3d   exp_rdc= %3d " % (len(rdcs), num_rdc)
        assert num_rdc == len(rdcs)
        sum_sq = 0.
        lns = []
        for i in range(num_rdc):
            diff =  exp_rdcs[i] - rdcs[i]
            sum_sq += diff * diff
            lns += [ "# %3d  resid= %3d  exp_rdc= %8.3f  cal_rdc= %8.3f  df= %8.3f " % (self.all_bp_ids[i],self.all_rdcs_resids[i], exp_rdcs[i], rdcs[i], exp_rdcs[i]-rdcs[i])]
        rmsd = sqrt(sum_sq/float(num_rdc))
        lns += [ "# Levfit: num_rdc= %3d  rmsd= %8.3f sum_sq= %8.3f " % (num_rdc, rmsd, sum_sq)  ]
        if vb: self.prnLns(lns)
        return rmsd, sum_sq,  lns  
    
    
    
    
    def setSimAngFilter(self):
        """  self.dup_angles : [(('1', '2'), [(-30.0, 30.0), (150.0, 210.0)]), (('1', '3'), [(-30.0, 30.0), (150.0, 210.0)]), (('2', '3'), [(-30.0, 30.0), (150.0, 210.0)])]"""
        if not hasattr(self, "dup_angles") or self.dup_angles is None:
            self.ANGFILT = 0
            return
        else:
            self.ANGFILT = 1
        self.dup_ang_filt = []
        for dup_ang in self.dup_angles:
            new_dp_angs = []
            for _dpa in dup_ang[1]:
                new_dp_angs.append((radians(_dpa[0]), radians(_dpa[1])))
            self.dup_ang_filt.append((dup_ang[0], new_dp_angs))
        #print self.dup_ang_filt
        
            
    def simFiltDupAng(self):
        for dp_aids, dp_angs in self.dup_ang_filt:
            dp1_id, dp2_id = dp_aids
            _ang = self.getVecAngle(self.all_dup_axis[dp1_id], self.all_dup_axis[dp2_id])
            if _ang < 0: _ang += pi2
            cnt = 0 
            for al, au in dp_angs:
                if _ang >= al and _ang <= au:
                    cnt = 1
                    break 
            if not cnt: return 0
        return 1 
            
        
    def simPM2DPR_DaR(self, parameters):
        num_dps = int((len(parameters) -2 )/ 4 ) #
        Da = parameters[-2]
        R = parameters[-1]
        dprs = []
        cnt = 0
        for i in range(num_dps):
            the = parameters[cnt]
            phi = parameters[cnt+1]
            rho = parameters[cnt+2]
            delta = parameters[cnt+3]
            dprs.append((the, phi, rho, delta))
            cnt += 4 
        return Da, R, dprs 
    
            
    def setSimDupAxle(self, parameters):
        """ return angles between two duplexs for duplex orientations """
        num_dps = int((len(parameters) -2 )/ 4 ) #
        Da = parameters[-2]
        R = parameters[-1]
        cnt = 0
        self.all_dup_axis = {}
        for i in range(num_dps):
            the = parameters[cnt]
            phi = parameters[cnt+1]
            cnt += 4
            dup_axle = self.thePhi2Axis(the, phi)
            self.all_dup_axis[self.duplex_ids[i]] = dup_axle
    
    
    def getSimDupAng(self):
        all_ang_dup = []
        num_dps = len(self.duplex_ids)
        for i in range(num_dps-1):
            dp1_id = self.duplex_ids[i]
            for j in range(i+1, num_dps):
                dp2_id = self.duplex_ids[j]
                _ang = self.getVecAngle(self.all_dup_axis[dp1_id], self.all_dup_axis[dp2_id])
                if _ang < 0: _ang += pi 
                all_ang_dup.append(degrees(_ang))
        #print all_ang_dup
        return tuple(all_ang_dup)
        
    
    
    def genRDC_Curve(self, pms, ndata):
        """ generate RDC curve with uniq slant angle """
        #pms = self.getTensor(helix=helix)
        #ndata = [ i*0.1 for i in range(5,120)]
        rdcs_uniq_slant = veglia(pms, None, ndata)
        lns = ["# No. n   rdc " ]
        for i in range(len(ndata)):
            lns += [ " %4d  %8.2f  %8.4f " %  (i, ndata[i], rdcs_uniq_slant[i])]
        self.prnLns(lns)
        #self.saveLines(lns, "arna_dp%s.dat" % helix, vb=1)
                        

        
class RDCDEG:
    """ RDC 4 degeneracy in duplex's referenec frame (The, Phi, Rho) """
    def prn4orients(self, o4, cmt=""):
        """ print RDC 4 orientations(The, Phi, Rho)"""
        print cmt 
        for i in range(len(o4)):
            _tpr = o4[i] #_tpr in o4:
            print "# i= %3d " % i, 
            print "# the= %6.1f phi= %6.1f rho= %6.1f" % _tpr 
            
    
    def gen4ThePhiRho(self, the, phi, rho, vb=0):
        """ generate 4 orientations by one set of the,phi, rho in units of deg """
        orientations = [self.__checkTPR(the, phi, rho),
                        self.__checkTPR(the, phi+180., rho),
                        self.__checkTPR(180.-the, 180-phi, rho + 180),
                        self.__checkTPR(180.-the, 360.-phi, rho + 180)]
        if vb: self.prn4orients(orientations)
        return orientations
    
    
    def gen4(self, the, phi, rho, vb=1, cmt=""):
        """ generate 4 orientations by one set of the,phi, rho in units of deg """
        o4 =  self.gen4ThePhiRho(the, phi, rho)
        if vb: self.prn4orients(o4, cmt=cmt)
        return o4 
        
    
    def __checkTPR(self, the, phi, rho):
        if phi < 0: phi += 360
        elif phi >360: phi -= 360
        if rho < 0: rho += 360
        elif rho> 360: rho -= 360
        return the, phi, rho
    
    def old_gen4(self, the, phi, rho):
        """ old way wrong generate 4 orientations by one set of the,phi, rho in units of deg """
        
        #if the > 90: the = 180. - the 
        #if phi > 180: phi = phi - 180. 
        #if rho > 180: 
        #    rho2 = rho - 180
        #else:
        #    rho2 = rho + 180 
        orientations = [(the, phi, rho),
                        (the, phi+180., rho),
                        (180.-the, 180-phi, rho + 180),
                        (180.-the, 360.-phi, rho + 180)]
        return orientations
