wavefunction.c

Go to the documentation of this file.

/*
** wavefunction.c
** 
** Made by (Johnny Q. Hacker)
** Login   <solarion@borkborkbork>
** 
** Stuff to deal with generalized wavefunctions.
Copyright (C) 2009 Joseph Pingenot

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Affero General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Affero General Public License for more details.

You should have received a copy of the GNU Affero General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
** Started on  Thu Mar 26 11:47:00 2009 Johnny Q. Hacker
** Last update Sun May 12 01:17:25 2002 Speed Blue
*/

#include <stdlib.h>
#include <math.h>
#include <libmodeling/wavefunction.h>
#include <libcommon/gcc_hints.h>

struct mstar_wavefunction_1d* mstar_wavefunction_1d_new(int N) {
  struct mstar_wavefunction_1d *wf = malloc(sizeof(struct mstar_wavefunction_1d));
  if(unlikely(wf == NULL)) return NULL;
  wf->points = malloc(N*sizeof(struct mstar_wavefunction_point_1d));
  if(unlikely(wf->points == NULL)) {
    free(wf);
    return NULL;
  }
  wf->N = N;
  return wf;
             }

void mstar_wavefunction_1d_free(struct mstar_wavefunction_1d* wf) {
  free(wf->points);
  free(wf);
}


struct mstar_wavefunction_1d* mstar_wavefunction_1d_new_gaussian(int N, double minx, double maxx, double sigma) {
  struct mstar_wavefunction_1d* wf = mstar_wavefunction_1d_new(N);
  if(unlikely(wf == NULL)) return NULL;
  double oosqrtsigma = 1.0/pow(sigma*M_PI, 0.25);
  double oosigma = 1.0/(sigma*2);
  for(int i=0; i<N; i++) {
    double x = (maxx-minx)/N*i + minx;
    (wf->points)[i].x = x;
    (wf->points)[i].re = oosqrtsigma*exp(-x*x*oosigma);
    (wf->points)[i].im = 0.0;
  }
  return wf;
}

struct mstar_wavefunction_1d* mstar_wavefunction_1d_from_LAPACK_normlized(const struct mstar_wavefunction_1d* wf) {
  struct mstar_wavefunction_1d* nwf = (struct mstar_wavefunction_1d*)malloc(sizeof(struct mstar_wavefunction_1d));
  if(nwf == NULL) return NULL;
  nwf->N = wf->N;
  nwf->points = (struct mstar_wavefunction_point_1d*)malloc(sizeof(struct mstar_wavefunction_point_1d)*nwf->N);
  if(nwf->points == NULL) {
    free(nwf);
    return NULL;
  }
  /*handle crazy cases gracefully*/
  if(unlikely(nwf->N < 3)) {
    if(nwf->N == 2) {
      nwf->points[0].x = wf->points[0].x;
      nwf->points[0].re = wf->points[0].re/sqrt(wf->points[1].x-wf->points[0].x);
      nwf->points[0].im = wf->points[0].im/sqrt(wf->points[1].x-wf->points[0].x);
      nwf->points[1].x = wf->points[1].x;
      nwf->points[1].re = wf->points[1].re/sqrt(wf->points[1].x-wf->points[0].x);
      nwf->points[1].im = wf->points[1].im/sqrt(wf->points[1].x-wf->points[0].x);
    }else{
      nwf->points[0].x = wf->points[0].x;
      nwf->points[0].re = wf->points[0].re;
      nwf->points[0].im = wf->points[0].im;
    }
  }else{
    /*Do average of forward and backward grid spaces, which necessitates special
     *  handling for the first and last points.*/
    nwf->points[0].x = wf->points[0].x;
    nwf->points[0].re = wf->points[0].re/sqrt(wf->points[1].x-wf->points[0].x);
    nwf->points[0].im = wf->points[0].im/sqrt(wf->points[1].x-wf->points[0].x);
    for(int i=1; i<(nwf->N-1); i++) {
      nwf->points[i].x = wf->points[i].x;
      nwf->points[i].re = wf->points[i].re/sqrt(0.5*(wf->points[i+1].x-wf->points[i-1].x));
      nwf->points[i].im = wf->points[i].im/sqrt(0.5*(wf->points[i+1].x-wf->points[i-1].x));
    }
    nwf->points[nwf->N-1].x = wf->points[nwf->N-1].x;
    nwf->points[nwf->N-1].re = wf->points[nwf->N-1].re/sqrt(wf->points[nwf->N-2].x-wf->points[nwf->N-1].x);
    nwf->points[nwf->N-1].im = wf->points[nwf->N-1].im/sqrt(wf->points[nwf->N-2].x-wf->points[nwf->N-1].x);
  }
  return nwf;
}

void mstar_wavefunction_1d_repack(const struct mstar_wavefunction_1d* wf, double **x, double **re, double **im, short int convert_lapack_to_density_normalization) {
  *x = malloc((wf->N)*sizeof(double));
  *re = malloc((wf->N)*sizeof(double));
  *im = malloc((wf->N)*sizeof(double));
  if(convert_lapack_to_density_normalization) {
    /*handle crazy cases gracefully*/
    if(unlikely(wf->N < 3)) {
      if(wf->N == 2) {
  (*x)[0] = wf->points[0].x;
  (*re)[0] = wf->points[0].re/sqrt(wf->points[1].x-wf->points[0].x);
  (*im)[0] = wf->points[0].im/sqrt(wf->points[1].x-wf->points[0].x);
  (*x)[1] = wf->points[1].x;
  (*re)[1] = wf->points[1].re/sqrt(wf->points[1].x-wf->points[0].x);
  (*im)[1] = wf->points[1].im/sqrt(wf->points[1].x-wf->points[0].x);
      }else{
  (*x)[0] = wf->points[0].x;
  (*re)[0] = wf->points[0].re;
  (*im)[0] = wf->points[0].im;
      }
    }else{
      /*Do average of forward and backward grid spaces, which necessitates special
       *  handling for the first and last points.*/
      (*x)[0] = wf->points[0].x;
      (*re)[0] = wf->points[0].re/sqrt(wf->points[1].x-wf->points[0].x);
      (*im)[0] = wf->points[0].im/sqrt(wf->points[1].x-wf->points[0].x);
      for(int i=1; i<(wf->N-1); i++) {
  (*x)[i] = wf->points[i].x;
  (*re)[i] = wf->points[i].re/sqrt(0.5*(wf->points[i+1].x-wf->points[i-1].x));
  (*im)[i] = wf->points[i].im/sqrt(0.5*(wf->points[i+1].x-wf->points[i-1].x));
      }
      (*x)[wf->N-1] = wf->points[wf->N-1].x;
      (*re)[wf->N-1] = wf->points[wf->N-1].re/sqrt(wf->points[wf->N-1].x-wf->points[wf->N-2].x);
      (*im)[wf->N-1] = wf->points[wf->N-1].im/sqrt(wf->points[wf->N-1].x-wf->points[wf->N-2].x);
    }
  }else{
    for(int i=0; i<wf->N; i++) {
      (*x)[i] = (wf->points)[i].x;
      (*re)[i] = (wf->points)[i].re;
      (*im)[i] = (wf->points)[i].im;
    }
  }
}