Message: Local Electric Field: GetFieldValue() problem Not Logged In (login)
 Next-in-Thread Next-in-Thread
 Next-in-Forum Next-in-Forum

Question Local Electric Field: GetFieldValue() problem 

Keywords: coordinate systems nonuniform electric field
Forum: Fields: Magnetic and Otherwise
Date: 15 Feb, 2008
From: Mario <Mario>

Hello Geant4 Gurus!!

I am using geant4.8.2.p01 currently and am implementing an electric field,
and later a magnetic field, to a specified local volume.  I have defined
my own electric field class, all of the code is below for that class.
The main concern is point[4] which comes in at the GetFieldValue() function.
I have displayed point[0], point[1], point[2] and calculated the magnitude
of this vector and in some instances it is larger than my entire world volume, 
which obviously makes no sense.  I guess I am curious about the coordinate 
system that point[4] is using... does it use the logicalVolume in which 
the fieldManager of the field belongs to as the coordinate system, or is it 
with respect to the global volume??  

Also, when I place this field_vol into the global and it overlaps with other
volumes already there or new ones that are placed on, it does not give me a
'overlap' warning, and the protons (which is what I am using) go through 
all volumes, overlapping or not, and it seems okay... is there something
special about a volume that is made of Vacuum and holds a field manager??

The main issue is what this point[4] is, because the EField is calculated
based on that number, so if I am using it incorrectly then there is a problem.
Thanks in advance for your help.

Cheers,
Mario 






HEADER:
#ifndef ChargeShell_hh
#define ChargeShell_hh 1

#include "G4Types.hh"
#include "G4ThreeVector.hh"
#include "G4ElectricField.hh"

class ChargeShell : public G4ElectricField
{
public: 

ChargeShell();
virtual ~ChargeShell() ;
virtual void GetFieldValue(const G4double point[4],
 G4double *field) const;

G4ThreeVector	CalcVolumeCharge(const G4double[4]) const;
G4ThreeVector	CalcSurfaceCharge(const G4double[4]) const;

G4double	GetLength()		{return 10.0*ro;};
G4ThreeVector	GetPos()		{return srcToShellPos;};

private:

	void ReadInputFile();

	G4double		ri;
	G4double		ro;
	G4double		Q;
	G4ThreeVector	srcToShellPos;

	G4double		fourPiEpsNot;
};
#endif


CC FILE:
#include "ChargeShell.hh"

extern int readFileVerbose;

ChargeShell::ChargeShell()
{
	ReadInputFile();
	fourPiEpsNot = 4.0*3.14159265359*8.85418781762e-12*farad/meter;
}

// ------------------------------------------------------------------------

ChargeShell::~ChargeShell()
{
}

void ChargeShell::GetFieldValue (const G4double point[4],G4double *fieldOut) const
{
	if(ri > ro){
		G4cout << "ri cannot be larger than ro" << G4endl;
		return;
	}
	
	G4ThreeVector fieldVal = G4ThreeVector(0.0,0.0,0.0);
	if(ro == ri){fieldVal = CalcSurfaceCharge(point);}
	else{fieldVal = CalcVolumeCharge(point);}

	G4cout << "fieldVal.x() = " << fieldVal.x()/(volt/m) << G4endl;
	G4cout << "fieldVal.y() = " << fieldVal.y()/(volt/m) << G4endl;
	G4cout << "fieldVal.z() = " << fieldVal.z()/(volt/m) << G4endl << G4endl;

	fieldOut[0]= 0.0;
	fieldOut[1]= 0.0;
	fieldOut[2]= 0.0;
	fieldOut[3]= fieldVal.x();
	fieldOut[4]= fieldVal.y();
	fieldOut[5]= fieldVal.z();
	return;
}

G4ThreeVector ChargeShell::CalcSurfaceCharge(const G4double p[4]) const 
{
	G4double x = p[0];
	G4double y = p[1];
	G4double z = p[2];
	G4double r = std::sqrt(x*x + y*y + z*z);

	G4double Er;
	G4cout << "r = " << r/mm << G4endl;
	if(r <= ro){Er = 0.0*volt/m;G4cout << "----->>>>>Er = 0<<<<<-----" << G4endl;}
	else{Er = Q/(fourPiEpsNot*r*r);}
	G4cout << "Er = " << Er/(volt/m) << G4endl;

	G4ThreeVector out = G4ThreeVector(Er*(x/r),Er*(y/r),Er*(z/r));
	return out;
}

G4ThreeVector ChargeShell::CalcVolumeCharge(const G4double p[4]) const
{
	G4double x = p[0];
	G4double y = p[1];
	G4double z = p[2];
	G4double r = std::sqrt(x*x + y*y + z*z);

	double riOverro = ri/ro;
	double riOverr = ri/r;

	G4double Er;
	if(r >= 0.0 && r <= ri){Er = 0.0*volt/m;}
	else if(r > ri && r < ro){
		Er = (Q/(fourPiEpsNot*ro*ro))*(r/ro)*(1-riOverr*riOverr*riOverr)/(1-riOverro*riOverro*riOverro);
	}
	else{Er = Q/(fourPiEpsNot*r*r);} 
	
	G4ThreeVector out = G4ThreeVector(Er*(x/r),Er*(y/r),Er*(z/r));
	return out;
}

void ChargeShell::ReadInputFile()
{
	//read in parameters from text file
	FILE *chargeShellInput;
	char* inputFile = "C:\\g4work\\PRSim\\src\\Param_Files\\Param_ChargeShellField.txt";
	chargeShellInput = fopen(inputFile,"rt");

	if(chargeShellInput == NULL){
		G4cout << "Param_ChargeShellField.txt did not come in correctly, check path" << G4endl;
		G4cout << "chargeShellInput = " << chargeShellInput << G4endl;
		G4cout << "errno = " << errno << G4endl;
		return;
	}	

	char readWholeLineBuff[500];
	float tmpFloat;
	char tmpString[50];

	//read line of text for ri
	fgets(readWholeLineBuff,500,chargeShellInput);
	if(readWholeLineBuff != NULL){
		fscanf(chargeShellInput,"%f\n",&tmpFloat);
		ri = tmpFloat*um;
		if(readFileVerbose > 0){
			G4cout << "readWholeLineBuff = " << readWholeLineBuff << G4endl;
			G4cout << ri/um << G4endl;
		}
	}
	else{G4cout << "Did not read in ri" << G4endl; return;}

	//read line of text for ro
	fgets(readWholeLineBuff,500,chargeShellInput);
	if(readWholeLineBuff != NULL){
		fscanf(chargeShellInput,"%f\n",&tmpFloat);
		ro = tmpFloat*um;
		if(readFileVerbose > 0){
			G4cout << "readWholeLineBuff = " << readWholeLineBuff << G4endl;
			G4cout << ro/um << G4endl;
		}
	}
	else{G4cout << "Did not read in ro" << G4endl; return;}

	//read line of text for Q
	fgets(readWholeLineBuff,500,chargeShellInput);
	fgets(readWholeLineBuff,500,chargeShellInput);
	if(readWholeLineBuff != NULL){
		fscanf(chargeShellInput,"%f\n",&tmpFloat);
		Q = tmpFloat*coulomb;
		if(readFileVerbose > 0){
			G4cout << "readWholeLineBuff = " << readWholeLineBuff << G4endl;
			G4cout << Q/coulomb << G4endl;
		}
	}
	else{G4cout << "Did not read in Q" << G4endl; return;}

	//read line of text for srcToShellPos
	fgets(readWholeLineBuff,500,chargeShellInput);
	if(readWholeLineBuff != NULL){
		fscanf(chargeShellInput,"%f",&tmpFloat);
		srcToShellPos.setX(tmpFloat*cm);
		fscanf(chargeShellInput,"%f",&tmpFloat);
		srcToShellPos.setY(tmpFloat*cm);
		fscanf(chargeShellInput,"%f\n",&tmpFloat);
		srcToShellPos.setZ(tmpFloat*cm);
		if(readFileVerbose > 0){
			G4cout << "readWholeLineBuff = " << readWholeLineBuff << G4endl;
			G4cout << srcToShellPos << G4endl;
		}
	}
	else{G4cout << "Did not read in srcToShellPos" << G4endl; return;}
	fclose(chargeShellInput);	
	return;
}

Inline Depth:
 1 1
 All All
Outline Depth:
 1 1
 2 2
 All All
Add message: (add)

1 Feedback: Re: Local Electric Field: GetFieldValue() problem   (Peter Gumplinger - 15 Feb, 2008)
(_ None: Re: Local Electric Field: GetFieldValue() problem   (Mario - 15 Feb, 2008)
(_ Feedback: Re: Local Electric Field: GetFieldValue() problem   (Peter Gumplinger - 18 Feb, 2008)
 Add Message Add Message
to: "Local Electric Field: GetFieldValue() problem"

 Subscribe Subscribe

This site runs SLAC HyperNews version 1.11-slac-98, derived from the original HyperNews


[ Geant 4 Home | Geant 4 HyperNews | Search | Request New Forum | Feedback ]