# -*- coding: utf-8 -*-
from sys import argv
from math import sin,cos,pi,sqrt,exp,cosh


minusfoursqrt3=-4.0*sqrt(3.0);sqrt3=sqrt(3.0);
sqrt3on2=sqrt(3.0)/2.0;twopionsqrt3=pi/sqrt3on2;twopion3=2.0*pi/3.0;fourpion3=2.0*twopion3

if len(argv)!=8:
  print "Give me a file, f, g, h, a, numpoints, and R0\n"
  exit(1);

def InitMatrix(argv,g,h,a):
  movingparts=open(argv,"r").readlines();movingparts.pop(0)
  MatrixVals={"qx":[],"qy":[],"D0":[],"D11":[],"D12":[]};

  for l in movingparts:
    lsplit=l.split()
    MatrixVals["qx"].append(float(lsplit[0]));MatrixVals["qy"].append(float(lsplit[1]))
    MatrixVals["D0"].append(g*float(lsplit[2])+h*float(lsplit[3]));
    MatrixVals["D11"].append(h*float(lsplit[4]));MatrixVals["D12"].append(h*float(lsplit[5]))
  
  pnum=len(MatrixVals["qx"])
  
  return MatrixVals

def CalculateMeanForce(MatrixVals,f,a,num,R0):
  dR=2.0*pi*(0.25-R0)/num;
  
  R=2.0*pi*R0;pnum=len(MatrixVals["qx"]);Prefac=2.0*(pi/a)**2;dqx=2.0*pi/sqrt(pnum);dqy=4.0*pi/sqrt(3.0*pnum);twoonsqrt3=2.0/sqrt(3.0);onethird=1.0/3.0;

  fp=open("MeanForce.dat","w");fp.write("#Rbar (1) Sigmax^2 (2) Sigmay^2 (3) Sigmaxy (4) Fsub (5)\n")
  for i in range(num+1):
    Sigmax=0.0;Sigmay=0.0;Sigmaxy=0.0;
    fpartD0=f*(2.0+cos(R));fpartD11=-4.0*f*(sin(R/2))**2
    
    for j in range(pnum):
      newD0=fpartD0+MatrixVals["D0"][j];newD11=fpartD11+MatrixVals["D11"][j]
      determinant=newD0*(newD0+newD11)-MatrixVals["D12"][j]**2
      
      if(MatrixVals["qx"][j]!=0.0 or MatrixVals["qy"][j]!=0.0):
        Sigmax =Sigmax+ newD0/(pnum*determinant);
        Sigmay =Sigmay+ (newD0+newD11)/(pnum*determinant);
        Sigmaxy=Sigmaxy-MatrixVals["D12"][j]/(pnum*determinant);
      
      else:Sigmay =Sigmay+ 1.0/(newD0*pnum);
        
    fp.write(str(R/(2.0*pi))+' '+str(Sigmax/a**2)+' '+str(Sigmay/a**2)+' '+str(Sigmaxy/a**2)+' '+str(sin(R)*cosh(twoonsqrt3*Prefac*Sigmaxy)*exp(-Prefac*(Sigmax+onethird*Sigmay) ) )+'\n')
    
    R=R+dR

[f,g,h,a,numpoints,R0]=[float(argv[2]),float(argv[3]),float(argv[4]),float(argv[5]),int(argv[6]),float(argv[7])]
CalculateMeanForce(InitMatrix(argv[1],g,h,a),f,a,numpoints,R0)