src/gmxlib/3dview.c

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