From 383c7d0abd2ccdbe73326afb7899ecefec65a0fe Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Thu, 20 Sep 2007 09:00:00 +0000 Subject: [PATCH] newly added ordered dipol field calculations, still preliminary --- dipol_fields/ordered_dipol_fields/LSCO.inp | 40 ++ dipol_fields/ordered_dipol_fields/Makefile | 28 + dipol_fields/ordered_dipol_fields/README | 5 + .../ordered_dipol_fields/dipolar_ordered.cpp | 633 ++++++++++++++++++ 4 files changed, 706 insertions(+) create mode 100644 dipol_fields/ordered_dipol_fields/LSCO.inp create mode 100644 dipol_fields/ordered_dipol_fields/Makefile create mode 100644 dipol_fields/ordered_dipol_fields/README create mode 100644 dipol_fields/ordered_dipol_fields/dipolar_ordered.cpp diff --git a/dipol_fields/ordered_dipol_fields/LSCO.inp b/dipol_fields/ordered_dipol_fields/LSCO.inp new file mode 100644 index 0000000..3a1127c --- /dev/null +++ b/dipol_fields/ordered_dipol_fields/LSCO.inp @@ -0,0 +1,40 @@ +%-------------------------------------------------------------------------------- +% +% LSCO.inp +% +% input file to calculate magnetic dipol field distribution for an ordered state +% +% data from Japanese Journal of Applied Physics, Part 1 (1988) 27, 1132-1137 +% +% muon site from H.-U. Suter, Physica B 326 (2003) 329 +%-------------------------------------------------------------------------------- + +%-------------------------------------------------------------------------------- +% magnetic lattice related input +%-------------------------------------------------------------------------------- + 4 % number of magnetic atoms in the magnetic unit cell + 2 % number of magnetic types in the magnetic unit cell + + 0.0000 0.0000 0.0000 1 % x, y, z, magnetic type + 0.5000 0.5000 0.0000 2 % x, y, z, magnetic type + 0.0000 0.5000 0.5000 1 % x, y, z, magnetic type + 0.5000 0.0000 0.5000 2 % x, y, z, magnetic type + + 5.3535 5.4003 13.1441 % lattice parameters in angstroms (magnetic cell) + 90.0000 90.0000 90.0000 % lattice angles in degrees + + 200.0000 % Lorentz sphere radius in angstroms + + 0.6800 0.0000 0.0000 % magnetic moment for type 1 + -0.6800 0.0000 0.0000 % magnetic moment for type 2 + +%-------------------------------------------------------------------------------- +% what to be calculated specifically follows now +%-------------------------------------------------------------------------------- + + point % tag showing that a single muon position to be calculated + 0.1478 0.0000 0.1846 % muon site in lattice coordinates + +%-------------------------------------------------------------------------------- +% end +%-------------------------------------------------------------------------------- diff --git a/dipol_fields/ordered_dipol_fields/Makefile b/dipol_fields/ordered_dipol_fields/Makefile new file mode 100644 index 0000000..1858b71 --- /dev/null +++ b/dipol_fields/ordered_dipol_fields/Makefile @@ -0,0 +1,28 @@ +#--------------------------------------------------------------------------- +# Makefile for the dipol field calculations +# +# Author: Andreas Suter, 2007-07-11 +# +# $Id$ +#--------------------------------------------------------------------------- + + +# compiler +CC = g++ +INCPATH = -I/usr/lib/qt-3.3/mkspecs/default -I$(QTDIR)/include +CFLAGS = -Wall -O2 -g +LDFLAGS = -Wl,-rpath,/usr/lib + +# libraries +LIBS = -L$(QTDIR)/lib -lqt-mt -lbsd -lm -lutil + +all: dipolar_ordered + +dipolar_ordered: dipolar_ordered.o + $(CC) $(LDFLAGS) -o dipolar_ordered dipolar_ordered.o $(LIBS) + +dipolar_ordered.o: dipolar_ordered.cpp + $(CC) $(CFLAGS) $(INCPATH) -c -o dipolar_ordered.o dipolar_ordered.cpp + +clean: + rm -rf *.o *~ diff --git a/dipol_fields/ordered_dipol_fields/README b/dipol_fields/ordered_dipol_fields/README new file mode 100644 index 0000000..bd69b58 --- /dev/null +++ b/dipol_fields/ordered_dipol_fields/README @@ -0,0 +1,5 @@ +Andreas Suter, 2007/09/20 + +This is still extremly experimental and not documented well, still it works. + +More to come ... diff --git a/dipol_fields/ordered_dipol_fields/dipolar_ordered.cpp b/dipol_fields/ordered_dipol_fields/dipolar_ordered.cpp new file mode 100644 index 0000000..c1e729d --- /dev/null +++ b/dipol_fields/ordered_dipol_fields/dipolar_ordered.cpp @@ -0,0 +1,633 @@ +#include +using namespace std; + +#include +#include +#include + +#include +#include +#include +#include + +#define DIPOLAR_SUCCESS 1 +#define DIPOLAR_FILE_NOT_FOUND -1 +#define DIPOLAR_FILE_OPEN_FAILED -2 +#define DIPOLAR_WRONG_INPUT_FILE_FORMAT -3 + +#define CALC_POINT 0 +#define CALC_PLANE 1 +#define CALC_SPHERE 2 + +// The half the number of unit cells used to calculate the dipolar field. +// This should be absorbed into the input file in the future. +#define DIPOLAR_NO_UNIT_CELLS 50 + +#define PI 3.141592653589793238462643383279502884197 + +// where is this number from? Check units and comment! +#define UNITS 9.2741 + +//----------------------------------------------------------------------- +/** + *

+ */ +typedef struct { + int line_no; + QString data; +} input_data; + +//----------------------------------------------------------------------- +/** + *

+ */ +typedef struct { + double pos[3]; + int type; +} atom_pos; + +//----------------------------------------------------------------------- +/** + *

+ */ +typedef struct { + double component[3]; +} mag_moment; + +//----------------------------------------------------------------------- +/** + *

+ */ +typedef struct { + int type; + double dipoleMatrix[3][3]; + double field[3]; +} dipole_type_info; + +//----------------------------------------------------------------------- +/** + *

+ * + *

calcTag: 0=single muon position, 1=dipole field in a plain, 2=dipole field on a sphere + */ +typedef struct { + int calcTag; /// tag showing what to be calculated + int noAtoms; + int noMagTypes; + QValueVector atom; + double latticeParam[3]; + double angle[3]; + double lorentzSphereRadius; + QValueVector magMoment; + double muPos[3]; + double planeFix[3]; /// vector fixing a plane in the crystal coordinate system + double planeNormal[3]; /// normal plane vector + QValueVector dipole; +} dipole_info; + +//----------------------------------------------------------------------- +/** + *

+ */ +void dipolar_ordered_syntax() +{ + cout << endl << "syntax: dipolar_ordered input_file output_file [debug]"; + cout << endl << " input_file: input file name"; + cout << endl << " The structure of the input file is:"; + cout << endl << " details to be written ..."; + cout << endl; + cout << endl << " WARNING: all the coordinates (including the magnetic moments) needed"; + cout << endl << " to be given in crystal structure coordinates. This is trivial"; + cout << endl << " for structures with lattice angles all 90°, but might lead to"; + cout << endl << " errors/confusions for lattice angles not equal to 90°"; + cout << endl << " a typical example would be:"; + cout << endl << " %--------------------------------------------------------------------------------"; + cout << endl << " %"; + cout << endl << " % LSCO.inp"; + cout << endl << " %"; + cout << endl << " % input file to calculate magnetic dipol field distribution for an ordered state"; + cout << endl << " %"; + cout << endl << " % data from Japanese Journal of Applied Physics, Part 1 (1988) 27, 1132-1137 "; + cout << endl << " % "; + cout << endl << " % muon site from H.-U. Suter, Physica B 326 (2003) 329"; + cout << endl << " %--------------------------------------------------------------------------------"; + cout << endl << ""; + cout << endl << " %--------------------------------------------------------------------------------"; + cout << endl << " % magnetic lattice related input"; + cout << endl << " %--------------------------------------------------------------------------------"; + cout << endl << " 4 % number of magnetic atoms in the magnetic unit cell"; + cout << endl << " 2 % number of magnetic types in the magnetic unit cell"; + cout << endl << ""; + cout << endl << " 0.0000 0.0000 0.0000 1 % x, y, z, magnetic type"; + cout << endl << " 0.5000 0.5000 0.0000 2 % x, y, z, magnetic type"; + cout << endl << " 0.0000 0.5000 0.5000 1 % x, y, z, magnetic type"; + cout << endl << " 0.5000 0.0000 0.5000 2 % x, y, z, magnetic type"; + cout << endl << ""; + cout << endl << " 5.3535 5.4003 13.1441 % lattice parameters in angstroms (magnetic cell)"; + cout << endl << " 90.0000 90.0000 90.0000 % lattice angles in degrees"; + cout << endl << ""; + cout << endl << " 200.0000 % Lorentz sphere radius in angstroms"; + cout << endl << ""; + cout << endl << " 0.6800 0.0000 0.0000 % magnetic moment for type 1"; + cout << endl << " -0.6800 0.0000 0.0000 % magnetic moment for type 2"; + cout << endl << ""; + cout << endl << " %--------------------------------------------------------------------------------"; + cout << endl << " % what to be calculated specifically follows now "; + cout << endl << " %--------------------------------------------------------------------------------"; + cout << endl << ""; + cout << endl << " point % tag showing that a single muon position to be calculated"; + cout << endl << " 0.1478 0.0000 0.1846 % muon site in lattice coordinates"; + cout << endl << ""; + cout << endl << " %--------------------------------------------------------------------------------"; + cout << endl << " % end"; + cout << endl << " %--------------------------------------------------------------------------------"; + cout << endl << " output_file: output file name"; + cout << endl << endl; +} + +//----------------------------------------------------------------------- +/** + *

+ */ +int dipolar_ordered_read_input_file(char *file_name, dipole_info *info) +{ + QFile file(file_name); + input_data line_data; + QValueVector data; + + bool error = false; + char line[256]; + QString qline; + int status; + int param_switch = 0; + int line_no = 0; + atom_pos atom; + mag_moment magMoment; + + // check if the file exists + if (!file.exists()) + return DIPOLAR_FILE_NOT_FOUND; + + // open file + if (!file.open(IO_ReadOnly)) + return DIPOLAR_FILE_OPEN_FAILED; + + // read all the data + int comment_pos; + while (file.readLine(&line[0], sizeof(line)) && !file.atEnd()) { + qline = line; + line_no++; + // ignore empty and comment lines + if (qline.stripWhiteSpace().isEmpty() || qline.startsWith("%")) + continue; + // strip comments within data line + comment_pos = qline.find("%"); + if (comment_pos>0) + qline = qline.remove(comment_pos-1, qline.length()-comment_pos); + qline = qline.lower(); + line_data.line_no = line_no; + line_data.data = qline; + data.append(line_data); + } + + file.close(); + + // analyze data + QValueVector::Iterator it; + for (it = data.begin(); it != data.end(); ++it) { + switch (param_switch) { + case 0: // no of magnetic atoms in the magnetic unit cell + status = sscanf(it->data.latin1(), "%d", &info->noAtoms); + if (status != 1) // likely to have wrong input file format + error = true; + param_switch++; + break; + case 1: // no of magnetic types in the magnetic unit cell + status = sscanf(it->data.latin1(), "%d", &info->noMagTypes); + if (status != 1) // likely to have wrong input file format + error = true; + param_switch++; + break; + case 2: // get all the atomic positions and the associated magnetic type + status = sscanf(it->data.latin1(), "%lf %lf %lf %d", + &atom.pos[0], &atom.pos[1], &atom.pos[2], &atom.type); + if (status != 4) { // likely to have wrong input file format + error = true; + } else { + info->atom.append(atom); + } + // check if all atom position are read + if (info->atom.count() == (unsigned int)info->noAtoms) { + param_switch++; + } + break; + case 3: // lattice parameters + status = sscanf(it->data.latin1(), "%lf %lf %lf", + &info->latticeParam[0], &info->latticeParam[1], &info->latticeParam[2]); + if (status != 3) // likely to have wrong input file format + error = true; + param_switch++; + break; + case 4: // lattice angles + status = sscanf(it->data.latin1(), "%lf %lf %lf", + &info->angle[0], &info->angle[1], &info->angle[2]); + if (status != 3) // likely to have wrong input file format + error = true; + param_switch++; + break; + case 5: // Lorentz sphere radius in (A) + status = sscanf(it->data.latin1(), "%lf", &info->lorentzSphereRadius); + if (status != 1) // likely to have wrong input file format + error = true; + param_switch++; + break; + case 6: // get all components of the magnetic types + status = sscanf(it->data.latin1(), "%lf %lf %lf", + &magMoment.component[0], &magMoment.component[1], &magMoment.component[2]); + if (status != 3) { // likely to have wrong input file format + error = true; + } else { + info->magMoment.append(magMoment); + } + // check if all magnetic moments are read + if (info->magMoment.count() == (unsigned int)info->noMagTypes) { + param_switch++; + } + break; + case 7: // tag what to be calculated + if (it->data.contains("point") > 0) { + info->calcTag = CALC_POINT; + } else if (it->data.contains("plane") > 0) { + info->calcTag = CALC_PLANE; + } else if (it->data.contains("sphere") > 0) { + info->calcTag = CALC_SPHERE; + } else { + error = true; + } + param_switch++; + break; + case 8: // depends on the calculation tag + switch (info->calcTag) { + case CALC_POINT: // get the muon position + status = sscanf(it->data.latin1(), "%lf %lf %lf", + &info->muPos[0], &info->muPos[1], &info->muPos[2]); + if (status != 3) + error = true; + break; + case CALC_PLANE: + status = sscanf(it->data.latin1(), "%lf %lf %lf", + &info->planeFix[0], &info->planeFix[1], &info->planeFix[2]); + if (status != 3) + error = true; + break; + case CALC_SPHERE: + break; + default: + break; + } + param_switch++; + break; + case 9: // depends on the calculation tag + switch (info->calcTag) { + case CALC_PLANE: + status = sscanf(it->data.latin1(), "%lf %lf %lf", + &info->planeNormal[0], &info->planeNormal[1], &info->planeNormal[2]); + if (status != 3) + error = true; + break; + default: + break; + } + param_switch++; + break; + default: + break; + } + + if (error) { + cout << endl << "ERROR in line " << it->line_no; + return DIPOLAR_WRONG_INPUT_FILE_FORMAT; + } + } + + return DIPOLAR_SUCCESS; +} + +//----------------------------------------------------------------------- +/** + *

+ */ +void dipolar_ordered_point_calc(dipole_info *info) +{ + double distance_vector[3]; + double direction_vector[3]; + double distance; ///< distance from an atom to the muon + double distance2; ///< squared distance form an atom to the muon + double idist3; ///< 1/distance^3 + double dcos[3]; ///< directional cosines + + // init dipole vector structure + dipole_type_info dipole; + for (int ii=0; ii<3; ii++) { + for (int jj=0; jj<3; jj++) { + dipole.dipoleMatrix[ii][jj] = 0.0; // init dipole tensor + } + dipole.field[ii] = 0.0; // init dipole field + } + for (int i=0; inoMagTypes; i++) { + dipole.type = i+1; // init magnetic type + info->dipole.append(dipole); // append to the vector + } + + // calculate the directional cosines + for (int i=0; i<3; i++) { + dcos[i] = cos(info->angle[i]*PI/180.0); + } + + // triple loop over the x, y, z unit cells + for (int uc_x = -DIPOLAR_NO_UNIT_CELLS; uc_x <= DIPOLAR_NO_UNIT_CELLS; uc_x++) { + for (int uc_y = -DIPOLAR_NO_UNIT_CELLS; uc_y <= DIPOLAR_NO_UNIT_CELLS; uc_y++) { + for (int uc_z = -DIPOLAR_NO_UNIT_CELLS; uc_z <= DIPOLAR_NO_UNIT_CELLS; uc_z++) { + // iterate all atoms in the unit cell + QValueVector::iterator it; + for (it = info->atom.begin(); it != info->atom.end(); ++it) { + + // calculate the distance vector of the magnetic atom to the muon + distance_vector[0] = info->latticeParam[0]*((it->pos[0]+uc_x) - info->muPos[0]); // "x" component + distance_vector[1] = info->latticeParam[1]*((it->pos[1]+uc_y) - info->muPos[1]); // "y" component + distance_vector[2] = info->latticeParam[2]*((it->pos[2]+uc_z) - info->muPos[2]); // "z" component + + // calculate the real distance between the atom and the muon + distance2 = distance_vector[0]*distance_vector[0]+ + distance_vector[1]*distance_vector[1]+ + distance_vector[2]*distance_vector[2]+ + distance_vector[0]*distance_vector[1]*dcos[0]+ + distance_vector[0]*distance_vector[2]*dcos[1]+ + distance_vector[1]*distance_vector[2]*dcos[2]; + distance = sqrt(distance2); + + // calculate the direction vector + direction_vector[0] = distance_vector[0]/distance; + direction_vector[1] = distance_vector[1]/distance; + direction_vector[2] = distance_vector[2]/distance; + + // check if it is within the Lorentz sphere, if yes calculate the dipole field tensor + if (info->lorentzSphereRadius - distance > 0.0) { + idist3 = 1.0/(distance*distance2); + for (int ii=0; ii<3; ii++) { + for (int jj=0; jj<3; jj++) { + if (ii==jj) { // diagonal elements + info->dipole.at(it->type-1).dipoleMatrix[ii][jj] += idist3*(3*direction_vector[ii]*direction_vector[ii]-1); + } else { // off diagonal elements + info->dipole.at(it->type-1).dipoleMatrix[ii][jj] += idist3*(3*direction_vector[ii]*direction_vector[jj]); + } + } + } + } + + } + } + } + } + + // calculate the field + for (int i=0; inoMagTypes; i++) { // calculation for each magnetic type in the lattice + for (int k=0; k<3; k++) { + for (int l=0; l<3; l++) { + info->dipole.at(i).field[k] += info->dipole.at(i).dipoleMatrix[k][l]*info->magMoment.at(i).component[l]; + } + } + } +} + +//----------------------------------------------------------------------- +/** + *

+ */ +void dipolar_ordered_plane_calc(dipole_info *info) +{ + cout << endl << "SORRY: not yet implemented ..." << endl; +} + +//----------------------------------------------------------------------- +/** + *

+ */ +void dipolar_ordered_sphere_calc(dipole_info *info) +{ + cout << endl << "SORRY: not yet implemented ..." << endl; +} + +//----------------------------------------------------------------------- +/** + *

+ */ +int dipolar_ordered_write_output_file(char *input_file_name, char *output_file_name, dipole_info *info) +{ + char str[128]; + // copy the input file to the output file + sprintf(str, "cp %s %s", input_file_name, output_file_name); + system(str); + + // open the new output file and append the results + QFile file(output_file_name); + + if (!file.open(IO_WriteOnly | IO_Append)) + return DIPOLAR_FILE_OPEN_FAILED; + + QTextStream fout(&file); + + double field_tot[3]; + for (int i=0; i<3; i++) + field_tot[i] = 0.0; + + fout << endl << "% --------- begin output data -----------------"; + fout << endl; + fout << endl << " ALL THE COMPONENTS ARE GIVEN AS COMPONENTS OF THE CRYSTAL, i.e."; + fout << endl << " vec_B = B_a vec_a + B_b vec_b + B_c vec_c, etc., where e.g. vec_a is the"; + fout << endl << " normalized vector along the a-axis of the crystal etc."; + fout << endl << " This needs to be kept in mind, especially for triklin, monoklin, ..., structures"; + fout << endl; + QValueVector::iterator it; + int i=1; + for (it = info->dipole.begin(); it != info->dipole.end(); ++it) { + fout << endl; + fout << endl << " THE DIPOLAR MATRIX FOR THE TYPE " << i++ << " in (kG/mu_B)"; + fout << endl; + fout << endl << " D | 1 2 3"; + fout << endl << " --|-----------------------------------"; + for (int j=0; j<3; j++) { + sprintf(str, " %d | %+2.4f %+2.4f %+2.4f ", + j+1, UNITS*it->dipoleMatrix[j][0], UNITS*it->dipoleMatrix[j][1], UNITS*it->dipoleMatrix[j][2]); + fout << endl << str; + field_tot[j] += UNITS*it->field[j]; + } + fout << endl; + fout << endl << " trace = " << UNITS*(it->dipoleMatrix[0][0]+it->dipoleMatrix[1][1]+it->dipoleMatrix[2][2]); + fout << endl; + fout << endl << " the dipolar magnetic field generated by the type " << it->type << " ions is: "; + fout << UNITS*sqrt(it->field[0]*it->field[0]+it->field[1]*it->field[1]+it->field[2]*it->field[2]) << " (kG)"; + fout << endl << " with the components: " << UNITS*it->field[0] << ", " << UNITS*it->field[1] << ", " << UNITS*it->field[2] << " (kG)"; + fout << endl; + fout << endl << " -----------------------"; + } + fout << endl; + fout << endl << " THE TOTAL DIPOLAR MAGNETIC FIELD IS "; + fout << sqrt(field_tot[0]*field_tot[0]+field_tot[1]*field_tot[1]+field_tot[2]*field_tot[2]) << " (kG)"; + fout << endl << " with the components: " << field_tot[0] << ", " << field_tot[1] << ", " << field_tot[2] << " (kG)"; + fout << endl; + fout << endl << "% --------- end output data -----------------"; + + file.close(); + + return DIPOLAR_SUCCESS; +} + +//----------------------------------------------------------------------- +/** + *

+ */ +int main(int argc, char *argv[]) +{ + int debug = 0; + dipole_info *info; + int status; + int i; + + // check syntax + if (argc < 3) { + dipolar_ordered_syntax(); + return 0; + } + + if (argc == 4) + sscanf(argv[3], "%d", &debug); + + // allocate memory for the info structure + info = new dipole_info; + + // init info structure + info->noAtoms = 0; + info->noMagTypes = 0; + info->latticeParam[0] = 0.0; + info->latticeParam[1] = 0.0; + info->latticeParam[2] = 0.0; + info->angle[0] = 0.0; + info->angle[1] = 0.0; + info->angle[2] = 0.0; + info->muPos[0] = 0.0; + info->muPos[1] = 0.0; + info->muPos[2] = 0.0; + info->lorentzSphereRadius = 0.0; + + // read input file and fill info structure + status = dipolar_ordered_read_input_file(argv[1], info); + if (status != DIPOLAR_SUCCESS) { + switch (status) { + case DIPOLAR_FILE_NOT_FOUND: + cout << endl << "ERROR: couldn't find input file " << argv[1] << endl; + break; + case DIPOLAR_FILE_OPEN_FAILED: + cout << endl << "ERROR: couldn't open input file " << argv[1] << endl; + break; + case DIPOLAR_WRONG_INPUT_FILE_FORMAT: + cout << endl << "ERROR: wrong input file format ..." << endl; + break; + default: + cout << endl << "unknown error ..." << endl; + break; + } + dipolar_ordered_syntax(); + return 0; + } + + if (debug) { + cout << endl << "-----------------------"; + cout << endl << "[debug] no of magnetic atoms in the unit cell = " << info->noAtoms; + cout << endl << "[debug] no of magnetic types in the unit cell = " << info->noMagTypes; + cout << endl << "[debug] list of the atom position (A) and its type"; + cout << endl << "-----------------------"; + i = 1; + for (QValueVector::iterator it = info->atom.begin(); it != info->atom.end(); ++it) { + cout << endl << "[debug] " << i++ << ": "; + cout << it->pos[0] << ", "; + cout << it->pos[1] << ", "; + cout << it->pos[2] << ", "; + cout << it->type; + } + cout << endl << "-----------------------"; + cout << endl << "[debug] lattice parameters (A) : "; + cout << info->latticeParam[0] << ", "; + cout << info->latticeParam[1] << ", "; + cout << info->latticeParam[2]; + cout << endl << "[debug] angles (°): "; + cout << info->angle[0] << ", "; + cout << info->angle[1] << ", "; + cout << info->angle[2]; + cout << endl << "[debug] Lorentz sphere radius (A): "; + cout << info->lorentzSphereRadius; + cout << endl << "[debug] magnetic moments of the various types"; + cout << endl << "-----------------------"; + i = 1; + for (QValueVector::iterator it = info->magMoment.begin(); it != info->magMoment.end(); ++it) { + cout << endl << "[debug] " << i++ << ": "; + cout << it->component[0] << ", "; + cout << it->component[1] << ", "; + cout << it->component[2] << ", "; + } + cout << endl << "-----------------------"; + if (info->calcTag == CALC_POINT) { + cout << endl << "[debug] calculate the dipolar field in a single point of the lattice."; + cout << endl << "[debug] muon position in lattice coordinates: "; + cout << info->muPos[0] << ", "; + cout << info->muPos[1] << ", "; + cout << info->muPos[2]; + } else if (info->calcTag == CALC_PLANE) { + cout << endl << "[debug] calculate the dipolar field in a plane of the lattice."; + cout << endl << "[debug] vector to the plane: "; + cout << info->planeFix[0] << ", "; + cout << info->planeFix[1] << ", "; + cout << info->planeFix[2]; + cout << endl << "[debug] plane normal vector: "; + cout << info->planeNormal[0] << ", "; + cout << info->planeNormal[1] << ", "; + cout << info->planeNormal[2]; + } else if (info->calcTag == CALC_SPHERE) { + cout << endl << "[debug] calculate the dipolar field on a sphere."; + } + cout << endl << "-----------------------"; + cout << endl << endl; + } + + // do the calculus + switch (info->calcTag) { + case CALC_POINT: + dipolar_ordered_point_calc(info); + break; + case CALC_PLANE: + dipolar_ordered_plane_calc(info); + break; + case CALC_SPHERE: + dipolar_ordered_sphere_calc(info); + break; + default: + cout << endl << "SEVERE ERROR: you never should get to here!?!?!?! Sorry ;-(" << endl; + break; + } + + // write the output file + status = dipolar_ordered_write_output_file(argv[1], argv[2], info); + if (status != DIPOLAR_SUCCESS) { + cout << endl << "ERROR while trying to write to the output file " << argv[2]; + cout << endl << endl; + } + + // clean up + if (info) { + info->atom.clear(); + delete info; + info = 0; + } + + return 1; +}