qpms/apps/hexlattice_ewald.c

   1 // c99 -o ew_altin -DALTIN -Wall -I ../.. -O2 -ggdb -DLATTICESUMS32 hexlattice_ewald.c ../translations.c ../ewald.c ../ewaldsf.c ../gaunt.c ../lattices2d.c -lgsl -lm -lblas
   2 // c99 -o ew_kin -DALTIN -DKSTDIN -Wall -I ../.. -O2 -ggdb -DLATTICESUMS32 hexlattice_ewald.c ../translations.c ../ewald.c ../ewaldsf.c ../gaunt.c ../lattices2d.c -lgsl -lm -lblas
   3 // c99 -o ew -Wall -I ../.. -O2 -ggdb -DLATTICESUMS32 hexlattice_ewald.c ../translations.c ../ewald.c ../ewaldsf.c ../gaunt.c ../lattices2d.c -lgsl -lm -lblas
   4 #include <stdio.h>
   5 #include <math.h>
   6 #include <string.h>
   7 #include <qpms/translations.h>
   8 #include <qpms/lattices.h>
   9 #include <gsl/gsl_const_mksa.h>
  10
  11
  12
  13 #define MAXOMEGACOUNT 1000
  14 #define MAXKCOUNT 200 // 200 // serves as klist default buffer size if KSTDIN is defined
  15 #define KMINCOEFF 0.783 //0.9783 // 0.783 // not used if KSTDIN defined
  16 #define KMAXCOEFF 1.217 //1.0217 // 1.217 // not used if KSTDIN defined
  17 #define KLAYERS 20
  18 #define RLAYERS 20
  19
  20 const double s3 = 1.732050807568877293527446341505872366942805253810380628055;
  21
  22 // IMPORTANT: lattice properties here
  23 const qpms_y_t lMax = 3;
  24 const double REFINDEX = 1.52;
  25 const double LATTICE_H = 576e-9;
  26 static const double SCUFF_OMEGAUNIT = 3e14;
  27 static const double hbar = GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR;
  28 static const double eV = GSL_CONST_MKSA_ELECTRON_CHARGE;
  29 static const double c0 = GSL_CONST_MKSA_SPEED_OF_LIGHT;
  30 static const TriangularLatticeOrientation rs_orientation = TRIANGULAR_VERTICAL;
  31
  32 int main (int argc, char **argv) {
  33   const double LATTICE_A = s3*LATTICE_H;
  34   const double INVLATTICE_A = 4*M_PI / s3 / LATTICE_A;
  35
  36   if (argc < 3) abort();
  37
  38   // char *kfile = argv[2]; // not used
  39   char *outfile = argv[2];
  40   char *errfile = NULL;
  41   if (argc > 3)
  42     errfile = argv[3];
  43
  44 #ifdef ALTIN // omega is provided on command line
  45   char *omegastr = argv[1];
  46   const double scuffomega = strtod(omegastr, NULL);
  47 #else
  48   char *omegafile = argv[1];
  49   double scuffomegas[MAXOMEGACOUNT];
  50   FILE *f = fopen(omegafile, "r");
  51   size_t omegacount = 0;
  52   while (fscanf(f, "%lf", scuffomegas + omegacount) == 1){
  53     assert(omegacount < MAXOMEGACOUNT);
  54     ++omegacount;
  55   }
  56   fclose(f);
  57 #endif
  58
  59   /*f = fopen(kfile, "r");
  60   int kcount = 100;
  61   while (fscanf(f, "%lf %lf", &(klist[kcount].x), &(klist[kcount].y)) == 2) {
  62     assert(kcount < MAXKCOUNT);
  63     ++kcount;
  64   }
  65   fclose(f);
  66   */
  67
  68 #ifdef KSTDIN
  69   size_t kcount = 0;
  70   size_t klist_capacity = MAXKCOUNT;
  71   cart2_t *klist = malloc(sizeof(cart2_t) * klist_capacity);
  72   while (scanf("%lf %lf", &(klist[kcount].x), &(klist[kcount].y)) == 2) {
  73     ++kcount;
  74     if(kcount >= klist_capacity) {
  75       klist_capacity *= 2;
  76       klist = realloc(klist, sizeof(cart2_t) * klist_capacity);
  77       if (klist == NULL) abort();
  78     }
  79   }
  80 #else
  81   cart2_t klist[MAXKCOUNT];
  82   int kcount = MAXKCOUNT;
  83   for (int i = 0; i <  kcount; ++i) { // TODO this should depend on orientation...
  84     klist[i].x = 0;
  85     klist[i].y = (4.* M_PI / 3. / LATTICE_A) * (KMINCOEFF + (KMAXCOEFF-KMINCOEFF)/kcount*i);
  86   }
  87 #endif
  88
  89   const double refindex = REFINDEX;
  90   const double h = LATTICE_H;
  91   const double a = h * s3;
  92   const double unitcell_area = s3*a*a/2.;
  93   const double rec_a = 4*M_PI/s3/a;
  94
  95   //fprintf(stderr, "%.16g\n", rec_a);
  96
  97   triangular_lattice_gen_t *Rlg = triangular_lattice_gen_init(a, rs_orientation, true, 0);
  98   triangular_lattice_gen_extend_to_steps(Rlg, RLAYERS);
  99   triangular_lattice_gen_t *Klg = triangular_lattice_gen_init(rec_a, reverseTriangularLatticeOrientation(rs_orientation), true, 0);
 100   triangular_lattice_gen_extend_to_steps(Klg, KLAYERS);
 101
 102   points2d_rordered_t Rselwo0 = points2d_rordered_annulus(&(Rlg->ps), 0, false, INFINITY, false);
 103
 104   const point2d *Rpoints = Rlg->ps.base;
 105   size_t nR = Rlg->ps.r_offsets[Rlg->ps.nrs];
 106   // TODO automatic skip of the zero point directly in translations.c
 107   const point2d *Rpoints_wo0 = Rselwo0.base + Rselwo0.r_offsets[0];
 108   size_t nR_wo0 = Rselwo0.r_offsets[Rselwo0.nrs] - Rselwo0.r_offsets[0];
 109   assert(nR - nR_wo0 == 1);
 110   const point2d *Kpoints = Klg->ps.base;
 111   size_t nK = Klg->ps.r_offsets[Klg->ps.nrs];
 112
 113   const point2d pshift0 = {0, 0};
 114   point2d pshiftAB = {0, 0}, pshiftBA = {0,0};
 115   if(rs_orientation == TRIANGULAR_VERTICAL) { // CHECKSIGN
 116     pshiftAB.x = h;///2;
 117     pshiftBA.x = -h;///2;
 118   } else { // CHECKSIGN
 119     pshiftAB.y = -h;///2;
 120     pshiftBA.y = h;///2;
 121   }
 122
 123   qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, QPMS_NORMALISATION_POWER_CS); // vai POWER_CS?
 124
 125   FILE *out = fopen(outfile, "w");
 126   FILE *err = NULL;
 127   if (errfile)
 128     err = fopen(errfile, "w");
 129
 130 #ifndef ALTIN
 131   for (size_t omegai = 0; omegai < omegacount; ++omegai) {
 132     const double scuffomega = scuffomegas[omegai];
 133 #else
 134   {
 135 #endif
 136     const double omega = scuffomega * SCUFF_OMEGAUNIT;
 137     const double EeV = omega * hbar / eV;
 138     const double k0_vac = omega / c0;
 139     const double k0_eff  = k0_vac * refindex;
 140     const double eta = 5.224/a; // FIXME quite arbitrary, but this one should work
 141
 142     // indices : destpart (A/B-particle), srcpart (A/B-particle), coeff type (A/B- type), desty, srcy
 143     complex double W[2][2][2][c->nelem][c->nelem];
 144     double Werr[2][2][2][c->nelem][c->nelem];
 145
 146     for (size_t ki = 0; ki < kcount; ++ki) {
 147       cart2_t beta = klist[ki];
 148       memset(W, 0, sizeof(W));
 149       if(err)
 150         memset(Werr, 0, sizeof(Werr));
 151
 152       const ptrdiff_t deststride = &(W[0][0][0][1][0]) - &(W[0][0][0][0][0]);
 153       const ptrdiff_t srcstride = &(W[0][0][0][0][1]) - &(W[0][0][0][0][0]);
 154       assert (srcstride == 1 && deststride == c->nelem);
 155
 156       // A<-A
 157       qpms_trans_calculator_get_AB_arrays_e32_both_points_and_shift(c,
 158           &(W[0][0][0][0][0]), err ? &(Werr[0][0][0][0][0]) : NULL, // Adest, Aerr,
 159           &(W[0][0][1][0][0]), err ? &(Werr[0][0][1][0][0]) : NULL, // Bdest, Berr,
 160           deststride, srcstride,
 161           eta, k0_eff, unitcell_area, nR_wo0, Rpoints_wo0, nK, Kpoints, beta, pshift0
 162           );
 163       // B<-B
 164       for(qpms_y_t desty = 0; desty < c->nelem; ++desty)
 165         for (qpms_y_t srcy = 0; srcy < c->nelem; ++srcy)
 166           for (int t = 0; t < 2; ++t) {
 167           W[1][1][t][desty][srcy] = W[0][0][t][desty][srcy];
 168           if (err)
 169             Werr[1][1][t][desty][srcy] = Werr[0][0][t][desty][srcy];
 170           }
 171
 172       // A<-B
 173       qpms_trans_calculator_get_AB_arrays_e32_both_points_and_shift(c,
 174           &(W[0][1][0][0][0]), err ? &(Werr[0][1][0][0][0]) : NULL, // Adest, Aerr,
 175           &(W[0][1][1][0][0]), err ? &(Werr[0][1][1][0][0]) : NULL, // Bdest, Berr,
 176           deststride, srcstride,
 177           eta, k0_eff, unitcell_area, nR, Rpoints, nK, Kpoints, beta, pshiftAB
 178           );
 179       // B<-A
 180       qpms_trans_calculator_get_AB_arrays_e32_both_points_and_shift(c,
 181           &(W[1][0][0][0][0]), err ? &(Werr[1][0][0][0][0]) : NULL, // Adest, Aerr,
 182           &(W[1][0][1][0][0]), err ? &(Werr[1][0][1][0][0]) : NULL, // Bdest, Berr,
 183           deststride, srcstride,
 184           eta, k0_eff, unitcell_area, nR, Rpoints, nK, Kpoints, beta, pshiftBA
 185           );
 186       // TODO CHECK B<-A vs. A<-B relation
 187
 188       fprintf(out, "%.16g\t%.16g\t%.16g\t%.16g\t%.16g\t",
 189           scuffomega, EeV, k0_eff, beta.x, beta.y);
 190       if(err) fprintf(err, "%.16g\t%.16g\t%16g\t%.16g\t%.16g\t",
 191           scuffomega, EeV, k0_eff, beta.x, beta.y);
 192       size_t totalelems = sizeof(W) / sizeof(complex double);
 193       for (size_t i = 0; i < totalelems; ++i) {
 194         complex double w = ((complex double *)W)[i];
 195         fprintf(out, "%.16g\t%.16g\t", creal(w), cimag(w));
 196         if (err)
 197           fprintf(err, "%.3g\t", ((double *)Werr)[i]);
 198       }
 199       fputc('\n', out);
 200       if(err) fputc('\n', err);
 201     }
 202   }
 203
 204   fclose(out);
 205   if(err) fclose(err);
 206
 207 #ifdef KSTDIN
 208   free(klist);
 209 #endif
 210
 211   qpms_trans_calculator_free(c);
 212
 213   triangular_lattice_gen_free(Klg);
 214   triangular_lattice_gen_free(Rlg);
 215 }