qpms/apps/2dlattice_ewald.c

   1 // c99 -o ew_gen_kin -Wall -I ../.. -I ../../amos/ -O2 -ggdb -DQPMS_VECTORS_NICE_TRANSFORMATIONS -DLATTICESUMS32 2dlattice_ewald.c ../translations.c ../ewald.c ../ewaldsf.c ../gaunt.c ../lattices2d.c ../latticegens.c ../bessel.c -lgsl -lm -lblas ../../amos/libamos.a -lgfortran ../error.c
   2 #include <stdio.h>
   3 #include <math.h>
   4 #include <string.h>
   5 #define LATTICESUMS32
   6 #include <qpms/translations.h>
   7 #include <qpms/lattices.h>
   8 #include <gsl/gsl_const_mksa.h>
   9
  10 // Command line order:
  11 // outfile b1.x b1.y b2.x b2.y lMax scuffomega refindex npart part0.x part0.y [part1.x part1.y [...]]
  12 //
  13 // Standard input (per line):
  14 // k.x k.y
  15 //
  16 // Output data format (line):
  17 //
  18
  19
  20 #define MAXKCOUNT 200 // 200 // serves as klist default buffer size
  21 //#define KMINCOEFF 0.783 //0.9783 // 0.783 // not used if KSTDIN defined
  22 //#define KMAXCOEFF 1.217 //1.0217 // 1.217 // not used if KSTDIN defined
  23 #define KLAYERS 20
  24 #define RLAYERS 20
  25
  26 const double s3 = 1.732050807568877293527446341505872366942805253810380628055;
  27
  28 // IMPORTANT: lattice properties here
  29 const qpms_y_t lMax = 3;
  30 const double REFINDEX = 1.52;
  31 const double LATTICE_H = 576e-9;
  32 static const double SCUFF_OMEGAUNIT = 3e14;
  33 static const double hbar = GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR;
  34 static const double eV = GSL_CONST_MKSA_ELECTRON_CHARGE;
  35 static const double c0 = GSL_CONST_MKSA_SPEED_OF_LIGHT;
  36
  37 int main (int argc, char **argv) {
  38   //const double LATTICE_A = s3*LATTICE_H;
  39   //const double INVLATTICE_A = 4*M_PI / s3 / LATTICE_A;
  40
  41   if (argc < 12) abort();
  42
  43   char *outfile = argv[1];
  44   char *errfile = NULL; // Filename for the error estimate output; NOT USED
  45   cart2_t b1 = {strtod(argv[2], NULL), strtod(argv[3], NULL)},
  46           b2 = {strtod(argv[4], NULL), strtod(argv[5], NULL)};
  47   const qpms_l_t lMax = strtol(argv[6], NULL, 10); assert(lMax>0);
  48   const double scuffomega = strtod(argv[7], NULL);
  49   const double refindex = strtod(argv[8], NULL);
  50   const int npart = strtol(argv[9], NULL, 10);  assert(npart>0);
  51   assert(argc >= 2*npart + 10);
  52   assert(npart > 0);
  53   cart2_t part_positions[npart];
  54   for(int p = 0; p < npart; ++p) {
  55     part_positions[p].x = strtod(argv[10+2*p], NULL);
  56     part_positions[p].y = strtod(argv[10+2*p+1], NULL);
  57   }
  58
  59 //#ifdef KSTDIN
  60   size_t kcount = 0;
  61   size_t klist_capacity = MAXKCOUNT;
  62   cart2_t *klist = malloc(sizeof(cart2_t) * klist_capacity);
  63   while (scanf("%lf %lf", &(klist[kcount].x), &(klist[kcount].y)) == 2) {
  64     ++kcount;
  65     if(kcount >= klist_capacity) {
  66       klist_capacity *= 2;
  67       klist = realloc(klist, sizeof(cart2_t) * klist_capacity);
  68       if (klist == NULL) abort();
  69     }
  70   }
  71 //#else
  72 #if 0
  73   cart2_t klist[MAXKCOUNT];
  74   int kcount = MAXKCOUNT;
  75   for (int i = 0; i <  kcount; ++i) { // TODO this should depend on orientation...
  76     klist[i].x = 0;
  77     klist[i].y = (4.* M_PI / 3. / LATTICE_A) * (KMINCOEFF + (KMAXCOEFF-KMINCOEFF)/kcount*i);
  78   }
  79 #endif
  80
  81   const double unitcell_area = l2d_unitcell_area(b1, b2);
  82   l2d_reduceBasis(b1, b2, &b1, &b2);
  83
  84   // TODO more clever way of determining the cutoff
  85   const double a = sqrt(unitcell_area); // N.B. different meaning than before
  86   const double maxR = 25 * a;
  87   const double maxK = 25 * 2*M_PI/a;
  88
  89   qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, QPMS_NORMALISATION_POWER_CS); // vai POWER_CS?
  90
  91   FILE *out = fopen(outfile, "w");
  92   FILE *err = NULL;
  93   if (errfile)
  94     err = fopen(errfile, "w");
  95
  96   {
  97     const double omega = scuffomega * SCUFF_OMEGAUNIT;
  98     const double EeV = omega * hbar / eV;
  99     const double k0_vac = omega / c0;
 100     const double k0_eff  = k0_vac * refindex;
 101     const double eta = 5.224/a; // FIXME quite arbitrary, but this one should work
 102
 103     // indices : destpart (A/B-particle), srcpart (A/B-particle), coeff type (A/B- type), desty, srcy
 104     complex double W[npart][npart][2][c->nelem][c->nelem];
 105     double Werr[npart][npart][npart][c->nelem][c->nelem];
 106
 107     for (size_t ki = 0; ki < kcount; ++ki) {
 108       cart2_t beta = klist[ki];
 109       memset(W, 0, sizeof(W));
 110       if(err)
 111         memset(Werr, 0, sizeof(Werr));
 112
 113       const ptrdiff_t deststride = &(W[0][0][0][1][0]) - &(W[0][0][0][0][0]);
 114       const ptrdiff_t srcstride = &(W[0][0][0][0][1]) - &(W[0][0][0][0][0]);
 115       assert (srcstride == 1 && deststride == c->nelem);
 116
 117       for (size_t ps = 0; ps < npart; ++ps) for (size_t pd = 0; pd < npart; ++pd)
 118         // TODO optimize (calculate only once for each particle shift; especially if pd == ps)
 119         qpms_trans_calculator_get_AB_arrays_e32(c,
 120             &(W[pd][ps][0][0][0]), err ? &(Werr[pd][ps][0][0][0]) : NULL, // Adest, Aerr,
 121             &(W[pd][ps][1][0][0]), err ? &(Werr[pd][ps][1][0][0]) : NULL, // Bdest, Berr,
 122             deststride, srcstride,
 123             eta, k0_eff, b1, b2,
 124             beta,
 125             cart2_substract(part_positions[pd], part_positions[ps]),  // CHECKSIGN
 126             maxR, maxK
 127             );
 128         // TODO CHECK B<-A vs. A<-B relation
 129
 130       fprintf(out, "%.16g\t%.16g\t%.16g\t%.16g\t%.16g\t",
 131           scuffomega, EeV, k0_eff, beta.x, beta.y);
 132       if(err) fprintf(err, "%.16g\t%.16g\t%16g\t%.16g\t%.16g\t",
 133           scuffomega, EeV, k0_eff, beta.x, beta.y);
 134       size_t totalelems = sizeof(W) / sizeof(complex double);
 135       for (size_t i = 0; i < totalelems; ++i) {
 136         complex double w = ((complex double *)W)[i];
 137         fprintf(out, "%.16g\t%.16g\t", creal(w), cimag(w));
 138         if (err)
 139           fprintf(err, "%.3g\t", ((double *)Werr)[i]);
 140       }
 141       fputc('\n', out);
 142       if(err) fputc('\n', err);
 143     }
 144   }
 145
 146   fclose(out);
 147   if(err) fclose(err);
 148
 149 //#ifdef KSTDIN
 150   free(klist);
 151 //#endif
 152
 153   qpms_trans_calculator_free(c);
 154
 155 }