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