The transverse current response function code for 3D free electron gas

//This code is to calculation the 3D free electron gas transverse 
//current response function  

#include 
#include 
#include 
#include 
#include 
#include 
using namespace std;

const double pi=3.141592653589793;

//temperature: we consider T=0;
//units: hbar=1;m=1;kf=1;
//chemical potential:
const double kf=1.0;

double Resuscep(double Qx,double omega);
double Imsuscep(double Qx,double omega);

int main()
{
	double omega,omegamax,domega;
	int iomega,Nomega;
	double Qx;
	ofstream Qx1,Qx3;


	//fix Qx=1, vary omega:
	Qx=1.0;
	omegamax=2*Qx+Qx*Qx+2.0;Nomega=1000;
        domega=omegamax/Nomega;
	Qx1.open("Qx1.txt");

	for(iomega=1;iomega<=2*Nomega+1;iomega++)
	{
		omega=-omegamax+(iomega-1)*domega;
		Qx1<2*Q-Q*Q && omega<2*Q+Q*Q)
		{
			Imsus=1/(2*Q)+pow(omega-Q*Q,2)/8/pow(Q,3)*(2*log(fabs(omega-Q*Q)/(2*Q))-1);
		}
		else if(omega<-2*Q+Q*Q && omega>-2*Q-Q*Q)
                {
                        Imsus=-1/(2*Q)-pow(-omega-Q*Q,2)/8/pow(Q,3)*(2*log(fabs(-omega-Q*Q)/(2*Q))-1);
                }
		else
		{
			Imsus=0;
		}
	}
	if(Q>2.0)
	{
		if(omega>-2*Q+Q*Q && omega<2*Q+Q*Q)
                {
                        Imsus=1/(2*Q)+pow(omega-Q*Q,2)/8/pow(Q,3)*(2*log(fabs(omega-Q*Q)/(2*Q))-1);
                }   
                else if(omega<2*Q-Q*Q && omega>-2*Q-Q*Q)
                {
                        Imsus=-1/(2*Q)-pow(-omega-Q*Q,2)/8/pow(Q,3)*(2*log(fabs(-omega-Q*Q)/(2*Q))-1);
                }
                else
                {
                        Imsus=0;
                }
	}

	return Imsus;
}

//The function to compute the real part of Lindhard susceptibility
//by the Kramers-Kronig relations
 
double Resuscep(double Q,double omega)
{
        double omegaprime,omegaprimemax,domegaprime;
        double Resus;
	double delta=10e-5;
	int Nomegaprime,iomegaprime;
	double eps=10e-10;

        omegaprimemax=Q*Q+2*Q+delta;
	Nomegaprime=1000;
	domegaprime=omegaprimemax/Nomegaprime;

	Resus=0.0;
	if(fabs(omega)>omegaprimemax)
	{
		for(iomegaprime=-Nomegaprime;iomegaprime<=Nomegaprime;iomegaprime++)
		{
			omegaprime=iomegaprime*omegaprimemax/Nomegaprime;
			Resus=Resus+1.0/pi*domegaprime*Imsuscep(Q,omegaprime)/(omegaprime-omega);
		}
	}
	if(fabs(omega)