src/gmxlib/3dview.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43 #include "sysstuff.h"
  44 #include "smalloc.h"
  45 #include "macros.h"
  46 #include "physics.h"
  47 #include "3dview.h"
  48 #include "pbc.h"
  49 #include "vec.h"
  50
  51 #define N 4
  52
  53 void m4_op(mat4 m,rvec x,vec4 v)
  54 {
  55   int i;
  56
  57   for(i=0; (i<N); i++)
  58     v[i]=m[XX][i]*x[XX]+m[YY][i]*x[YY]+m[ZZ][i]*x[ZZ]+m[WW][i];
  59 }
  60
  61 void unity_m4(mat4 m)
  62 {
  63   int i,j;
  64
  65   for(i=0; (i<N); i++)
  66     for(j=0; (j<N); j++)
  67       if (i==j)
  68         m[i][j]=1.0;
  69       else
  70         m[i][j]=0.0;
  71 }
  72
  73 void print_m4(FILE *fp,const char *s,mat4 A)
  74 {
  75   int i,j;
  76
  77   if (fp) {
  78     fprintf(fp,"%s: ",s);
  79     for (i=0; i<N; i++) {
  80       fprintf(fp,"\t");
  81       for (j=0; j<N; j++)
  82         fprintf(fp,"%10.5f",A[i][j]);
  83       fprintf(fp,"\n");
  84     }
  85   }
  86 }
  87
  88 void print_v4(FILE *fp,char *s,int dim,real *a)
  89 {
  90   int j;
  91
  92   if (fp) {
  93     fprintf(fp,"%s: ",s);
  94     for (j=0; j<dim; j++)
  95       fprintf(fp,"%10.5f",a[j]);
  96     fprintf(fp,"\n");
  97   }
  98 }
  99
 100 void mult_matrix(mat4 A, mat4 B, mat4 C)
 101 {
 102   int i,j,k;
 103
 104   for (i=0; i<N; i++)
 105     for (j=0; j<N; j++) {
 106       A[i][j]=0;
 107       for(k=0; (k<N); k++)
 108         A[i][j]+=B[i][k]*C[k][j];
 109     }
 110 }
 111
 112 void rotate(int axis, real angle, mat4 A)
 113 {
 114   unity_m4(A);
 115
 116   switch (axis) {
 117   case XX:
 118     A[YY][YY] =  cos(angle);
 119     A[YY][ZZ] = -sin(angle);
 120     A[ZZ][YY] =  sin(angle);
 121     A[ZZ][ZZ] =  cos(angle);
 122     break;
 123   case YY:
 124     A[XX][XX] =  cos(angle);
 125     A[XX][ZZ] =  sin(angle);
 126     A[ZZ][XX] = -sin(angle);
 127     A[ZZ][ZZ] =  cos(angle);
 128     break;
 129   case ZZ:
 130     A[XX][XX] =  cos(angle);
 131     A[XX][YY] = -sin(angle);
 132     A[YY][XX] =  sin(angle);
 133     A[YY][YY] =  cos(angle);
 134     break;
 135   default:
 136     gmx_fatal(FARGS,"Error: invalid axis: %d",axis);
 137   }
 138 }
 139
 140 void translate(real tx, real ty, real tz, mat4 A)
 141 {
 142   unity_m4(A);
 143   A[3][XX] = tx;
 144   A[3][YY] = ty;
 145   A[3][ZZ] = tz;
 146 }
 147
 148 static void set_scale(t_3dview *view,real sx, real sy)
 149 {
 150   view->sc_x=sx;
 151   view->sc_y=sy;
 152 }
 153
 154 void calculate_view(t_3dview *view)
 155 {
 156 #define SMALL 1e-6
 157   mat4 To,Te,T1,T2,T3,T4,T5,N1,D1,D2,D3,D4,D5;
 158   real dx,dy,dz,l,r;
 159
 160   /* eye center */
 161   dx=view->eye[XX];
 162   dy=view->eye[YY];
 163   dz=view->eye[ZZ];
 164   l = sqrt(dx*dx+dy*dy+dz*dz);
 165   r = sqrt(dx*dx+dy*dy);
 166 #ifdef DEBUG
 167   print_v4(debug,"eye",N,view->eye);
 168   printf("del: %10.5f%10.5f%10.5f l: %10.5f, r: %10.5f\n",dx,dy,dz,l,r);
 169 #endif
 170   if (l < SMALL)
 171     gmx_fatal(FARGS,"Error: Zero Length Vector - No View Specified");
 172   translate((real)(-view->origin[XX]),
 173             (real)(-view->origin[YY]),(real)(-view->origin[ZZ]),To);
 174   translate((real)(-view->eye[XX]),
 175             (real)(-view->eye[YY]),(real)(-view->eye[ZZ]),Te);
 176
 177   unity_m4(T2);
 178   T2[YY][YY]=0, T2[YY][ZZ]=-1, T2[ZZ][YY]=1, T2[ZZ][ZZ]=0;
 179
 180   unity_m4(T3);
 181   if (r > 0)
 182     T3[XX][XX]=-dy/r, T3[XX][ZZ]=dx/r, T3[ZZ][XX]=-dx/r, T3[ZZ][ZZ]=-dy/r;
 183
 184   unity_m4(T4);
 185   T4[YY][YY]=r/l, T4[YY][ZZ]=dz/l, T4[ZZ][YY]=-dz/l, T4[ZZ][ZZ]=r/l;
 186
 187   unity_m4(T5);
 188   T5[ZZ][ZZ]=-1;
 189
 190   unity_m4(N1);
 191   /* N1[XX][XX]=4,N1[YY][YY]=4; */
 192
 193   mult_matrix(T1,To,view->Rot);
 194   mult_matrix(D1,Te,T2);
 195   mult_matrix(D2,T3,T4);
 196   mult_matrix(D3,T5,N1);
 197   mult_matrix(D4,T1,D1);
 198   mult_matrix(D5,D2,D3);
 199
 200   mult_matrix(view->proj,D4,D5);
 201
 202 #ifdef DEBUG
 203   print_m4(debug,"T1",T1);
 204   print_m4(debug,"T2",T2);
 205   print_m4(debug,"T3",T3);
 206   print_m4(debug,"T4",T4);
 207   print_m4(debug,"T5",T5);
 208   print_m4(debug,"N1",N1);
 209   print_m4(debug,"Rot",view->Rot);
 210   print_m4(debug,"Proj",view->proj);
 211 #endif
 212 }
 213
 214 gmx_bool zoom_3d(t_3dview *view,real fac)
 215 {
 216   real dr;
 217   real bm,dr1,dr2;
 218   int  i;
 219
 220   dr2=0;
 221   for(i=0; (i<DIM); i++) {
 222     dr=view->eye[i];
 223     dr2+=dr*dr;
 224   }
 225   dr1=sqrt(dr2);
 226   if (fac < 1) {
 227     bm=max(norm(view->box[XX]),max(norm(view->box[YY]),norm(view->box[ZZ])));
 228     if (dr1*fac < 1.1*bm) /* Don't come to close */
 229       return FALSE;
 230   }
 231
 232   for(i=0; (i<DIM); i++)
 233     view->eye[i]*=fac;
 234   calculate_view(view);
 235   return TRUE;
 236 }
 237
 238 void init_rotate_3d(t_3dview *view)
 239 {
 240   real rot=DEG2RAD*15;
 241   int i;
 242
 243   for(i=0; (i<DIM); i++) {
 244     rotate(i,        rot ,view->RotP[i]);
 245     rotate(i,(real)(-rot),view->RotM[i]);
 246 #ifdef DEBUG
 247     print_m4(debug,"RotP",view->RotP[i]);
 248     print_m4(debug,"RotM",view->RotM[i]);
 249 #endif
 250   }
 251 }
 252
 253
 254 void rotate_3d(t_3dview *view,int axis,gmx_bool bPositive)
 255 {
 256   int  i,j;
 257   mat4 m4;
 258
 259   if (bPositive)
 260     mult_matrix(m4,view->Rot,view->RotP[axis]);
 261   else
 262     mult_matrix(m4,view->Rot,view->RotM[axis]);
 263   for(i=0; (i<N); i++)
 264     for(j=0; (j<N); j++)
 265     view->Rot[i][j]=m4[i][j];
 266
 267   calculate_view(view);
 268 }
 269
 270 void translate_view(t_3dview *view,int axis,gmx_bool bPositive)
 271 {
 272 #ifdef DEBUG
 273   printf("Translate called\n");
 274 #endif
 275   if (bPositive)
 276     view->origin[axis]+=view->box[axis][axis]/8;
 277   else
 278     view->origin[axis]-=view->box[axis][axis]/8;
 279   calculate_view(view);
 280 }
 281
 282 void reset_view(t_3dview *view)
 283 {
 284   int  i;
 285
 286 #ifdef DEBUG
 287   printf("Reset view called\n");
 288 #endif
 289   set_scale(view,4.0,4.0);
 290   clear_rvec(view->eye);
 291   calc_box_center(view->ecenter,view->box,view->origin);
 292   view->eye[ZZ]=3.0*max(view->box[XX][XX],view->box[YY][YY]);
 293   zoom_3d(view,1.0);
 294   view->eye[WW]=view->origin[WW]=0.0;
 295
 296   /* Initiate the matrix */
 297   unity_m4(view->Rot);
 298   calculate_view(view);
 299
 300   init_rotate_3d(view);
 301 }
 302
 303 t_3dview *init_view(matrix box)
 304 {
 305   t_3dview *view;
 306   int      i,j;
 307
 308   snew(view,1);
 309
 310   /* Copy parameters into variables */
 311   for(i=0; (i<DIM); i++)
 312     for(j=0; (j<DIM); j++)
 313       view->box[i][j]=box[i][j];
 314
 315   view->ecenter = ecenterDEF;
 316
 317   reset_view(view);
 318
 319   return view;
 320 }
 321