1 /* Copyright (C) 2003-2007 CAMP
2 * Copyright (C) 2007-2008 CAMd
3 * Please see the accompanying LICENSE file for further information. */
4
5 #include <Python.h>
6 #define PY_ARRAY_UNIQUE_SYMBOL GPAW_ARRAY_API
7 #define NO_IMPORT_ARRAY
8 #include <numpy/arrayobject.h>
9 #include "extensions.h"
10
11 int write_plt_file(char *fname,
12 int nx, int ny, int nz,
13 double x0, double y0, double z0,
14 double dx, double dy, double dz,
15 double *grid);
16
17 /* write grid to binary plt (gOpenMol) plot file */
WritePLT(PyObject * self,PyObject * args)18 PyObject* WritePLT(PyObject *self, PyObject *args)
19 {
20 char* fname; /* file name */
21 PyArrayObject* ho; /* grid spacings */
22 PyArrayObject* go; /* grid to write */
23 if (!PyArg_ParseTuple(args, "sOO", &fname, &ho, &go))
24 return NULL;
25
26 /* must be 3D */
27 if(PyArray_NDIM(go) != 3) return NULL;
28
29 double* g = DOUBLEP(go);
30 double* h = DOUBLEP(ho);
31
32 write_plt_file(fname,
33 PyArray_DIM(go, 0),
34 PyArray_DIM(go, 1),
35 PyArray_DIM(go, 2),
36 0.,0.,0.,
37 h[0],h[1],h[2],
38 g);
39 Py_RETURN_NONE;
40 }
41
42 /* -----------------------------------------------------------------
43 * write grid to binary plt (gOpenMol) plot file
44 *
45 * x0, dx etc are assumed to be atomic units
46 * the grid is assumed to be in the format:
47 * grid(ix,iy,iz) = grid[ ix + ( iy + iz*ny )*nx ];
48 * where ix=0..nx-1 etc
49 */
50
51 /* stolen from pltfile.c */
52 #define FWRITE(value , size) { \
53 Items = fwrite(&value, size , 1L , Output_p);\
54 if(Items < 1) {\
55 printf("?ERROR - in writing contour file (*)\n");\
56 return(1);}}
57
write_plt_file(char * fname,int nx,int ny,int nz,double x0,double y0,double z0,double dx,double dy,double dz,double * grid)58 int write_plt_file(char *fname,
59 int nx, int ny, int nz,
60 double x0, double y0, double z0,
61 double dx, double dy, double dz,
62 double *grid) {
63 FILE *Output_p;
64 static int Items;
65 float scale,zmin,zmax,ymin,ymax,xmin,xmax,val;
66 int rank,TypeOfSurface;
67 int ix,iy,iz,indx;
68 double norm,sum,dV;
69
70 Output_p = fopen(fname,"wb");
71
72 /* see http://www.csc.fi/gopenmol/developers/plt_format.phtml */
73
74 #define au_A 0.52917725
75 scale = au_A; /* atomic length in Angstroem */
76
77 rank=3; /* always 3 */
78 FWRITE(rank , sizeof(int));
79 TypeOfSurface=4; /* arbitrary */
80 FWRITE(TypeOfSurface , sizeof(int));
81 FWRITE(nz , sizeof(int));
82 FWRITE(ny , sizeof(int));
83 FWRITE(nx , sizeof(int));
84 zmin= scale * ((float) z0);
85 zmax= scale * ((float) z0+(nz-1)*dz);
86 /* float zmax=(float) z0+nz*dz; */
87 FWRITE(zmin , sizeof(float));
88 FWRITE(zmax , sizeof(float));
89 ymin= scale * ((float) y0);
90 ymax= scale * ((float) y0+(ny-1)*dy);
91 /* float ymax=(float) y0+ny*dy; */
92 FWRITE(ymin , sizeof(float));
93 FWRITE(ymax , sizeof(float));
94 xmin= scale * ((float) x0);
95 xmax= scale * ((float) x0+(nx-1)*dx);
96 /* float xmax=(float) x0+nx*dx; */
97 FWRITE(xmin , sizeof(float));
98 FWRITE(xmax , sizeof(float));
99
100 indx=0;
101 norm = 0;
102 sum=0;
103 dV=dx*dy*dz;
104 for(iz=0;iz<nz;iz++)
105 for(iy=0;iy<ny;iy++)
106 for(ix=0;ix<nx;ix++) {
107 val = (float) grid[indx];
108 sum += val;
109 norm += val*val;
110 FWRITE(val , sizeof(float));
111 indx++;
112 }
113
114 fclose(Output_p);
115
116 printf("#<write_plt_file> %s written (sum=%g,norm=%g)\n",
117 fname,sum*dV,norm*dV);
118
119 return 0;
120 }
121
122