Ignore:
Timestamp:
Jun 21, 2008, 8:18:51 PM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
bd3839
Parents:
1b192d
Message:

Finally, NanoCreator is working correctly!

Sheet has correct lengths in BOTH directions, hence tube is nice and alignment works, and hence torus works and alignment works also!

File:
1 edited

Legend:

Unmodified
Added
Removed

  • util/src/NanoCreator.c

    r1b192d re2f3fde  
    55#include <time.h>
    66
     7#include <gsl/gsl_matrix.h>
     8#include <gsl/gsl_vector.h>
     9#include <gsl/gsl_eigen.h>
     10#include <gsl/gsl_blas.h>
     11
    712#include "NanoCreator.h"
     13
    814
    915// ---------------------------------- F U N C T I O N S ----------------------------------------------
    … …  
    879885        Tubevector[axis[2]][i] = Vector[axis[2]][i];
    880886      }
    881 //      tmp =  Norm(Tubevector[axis[1]]);
    882 //      for (i=0;i<NDIM;i++)
    883 //        Tubevector[axis[1]][i] *= Norm(Vector[axis[1]])*sqrt(3)*sqrt(chiral[0]*chiral[0]+chiral[1]*chiral[1]+chiral[0]*chiral[1])/(double)ggT/tmp;
    884       fprintf(stdout, "Orthogonal Tubevector axis[0->1] %lg\n", Projection(Tubevector[axis[0]], Tubevector[axis[1]]));
    885       fprintf(stdout, "Orthogonal Tubevector axis[0->2] %lg\n", Projection(Tubevector[axis[0]], Tubevector[axis[2]]));
    886       fprintf(stdout, "Orthogonal Tubevector axis[1->2] %lg\n", Projection(Tubevector[axis[1]], Tubevector[axis[2]]));
     887      // here we assume, that Vector[axis[2]] is along z direction!
     888      gsl_matrix *M = gsl_matrix_alloc(2,2);
     889      gsl_matrix *C = gsl_matrix_alloc(2,2);
     890      gsl_matrix *evec = gsl_matrix_alloc(2,2);
     891      gsl_vector *eval = gsl_vector_alloc(2);
     892      gsl_vector *v = gsl_vector_alloc(2);
     893      gsl_vector *u = gsl_vector_alloc(2);
     894      gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(2);
     895      gsl_matrix_set(C, 0,0, Vector[axis[0]][0]);
     896      gsl_matrix_set(C, 1,0, Vector[axis[0]][1]);
     897      gsl_matrix_set(C, 0,1, Vector[axis[1]][0]);
     898      gsl_matrix_set(C, 1,1, Vector[axis[1]][1]);
     899      gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0, C, C, 0.0, M);
     900      fprintf(stdout, "M: \t%lg\t%lg\n\t%lg\t%lg\n", gsl_matrix_get(M,0,0), gsl_matrix_get(M,0,1), gsl_matrix_get(M,1,0), gsl_matrix_get(M,1,1));
     901      gsl_eigen_symmv(M, eval, evec, w);
     902      gsl_eigen_symmv_sort(eval,evec,GSL_EIGEN_SORT_ABS_DESC);
     903      fprintf(stdout, "Eigenvalues: %lg\t%lg\n", gsl_vector_get(eval,0), gsl_vector_get(eval,1));
     904      fprintf(stdout, "Eigenvectors: \t%lg\t%lg\n\t\t%lg\t%lg\n", gsl_matrix_get(evec,0,0), gsl_matrix_get(evec,0,1), gsl_matrix_get(evec,1,0), gsl_matrix_get(evec,1,1));
     905      gsl_matrix_set(M, 0,0, 0.);
     906      gsl_matrix_set(M, 1,0, 1.);
     907      gsl_matrix_set(M, 0,1, -gsl_vector_get(eval,1)/gsl_vector_get(eval,0));
     908      gsl_matrix_set(M, 1,1, 0.);
     909      gsl_vector_set(v,0,(double)chiral[0]);
     910      gsl_vector_set(v,1,(double)chiral[1]);
     911      gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0, evec, M, 0.0, C);
     912      gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0, C, evec, 0.0, M);
     913      fprintf(stdout, "M: \t%lg\t%lg\n\t%lg\t%lg\n", gsl_matrix_get(M,0,0), gsl_matrix_get(M,0,1), gsl_matrix_get(M,1,0), gsl_matrix_get(M,1,1));
     914      gsl_blas_dgemv(CblasNoTrans, 1.0, M, v, 0.0, u);
     915      fprintf(stdout, "Looking for factor to integer...\n");
     916      for(i=1;i<(chiral[0]+chiral[1])*(chiral[0]+chiral[1]);i++) {
     917        x1 = gsl_vector_get(u,0)*(double)i;
     918        x2 = gsl_vector_get(u,1)*(double)i;
     919        x3 =
     920        fprintf(stdout, "%d: %d\t%d vs. %lg\t%lg\n",i, ((int)(x1+x1/fabs(x1)*.5)), ((int)(x2+x2/fabs(x2)*.5)), (x1), (x2));
     921        if (( fabs( ((int)(x1+x1/fabs(x1)*.5)) - (x1) ) < 1e-3) && ( fabs( ((int)(x2+x2/fabs(x2)*.5)) - (x2) ) < 1e-3 )) {
     922          gsl_blas_dscal((double)i, u);
     923          break;
     924        }
     925      }
     926      fprintf(stdout, "(c,d) = (%lg,%lg)\n",gsl_vector_get(u,0), gsl_vector_get(u,1));
     927
     928      // get length
     929      double x[NDIM];
     930      for (i=0;i<NDIM;i++)
     931        x[i] = gsl_vector_get(u,0) * Vector[axis[0]][i] + gsl_vector_get(u,1) * Vector[axis[1]][i];
     932      angle = Norm(x)/Norm(Tubevector[axis[1]]) ;//ScalarProduct(x,Tubevector[axis[1]])/Norm(Tubevector[axis[1]]);
     933      fprintf(stdout, "angle is %lg\n", angle);
     934      for (i=0;i<NDIM;i++) {
     935        Tubevector[axis[1]][i] = gsl_vector_get(u,0) * Vector[axis[0]][i] + gsl_vector_get(u,1) * Vector[axis[1]][i];
     936      }
     937
     938      // Probe
     939      gsl_matrix_set(M, 0,0, Vector[axis[0]][0]);
     940      gsl_matrix_set(M, 1,0, Vector[axis[0]][1]);
     941      gsl_matrix_set(M, 0,1, Vector[axis[1]][0]);
     942      gsl_matrix_set(M, 1,1, Vector[axis[1]][1]);
     943      gsl_vector_set(v,0,(double)chiral[0]);
     944      gsl_vector_set(v,1,(double)chiral[1]);
     945      gsl_blas_dgemv(CblasNoTrans, 1.0, M, u, 0.0, eval);
     946      gsl_blas_dgemv(CblasTrans, 1.0, M, v, 0.0, u);
     947      x1=1.;
     948      gsl_blas_ddot(u,eval,&x1);
     949      fprintf(stdout, "Testing (c,d): (a,b) M^t M (c,d)^t = 0 ? : %lg\n", x1);
     950
     951      gsl_matrix_free(M);
     952      gsl_matrix_free(C);
     953      gsl_matrix_free(evec);
     954      gsl_vector_free(eval);
     955      gsl_vector_free(v);
     956      gsl_vector_free(u);
     957      gsl_eigen_symmv_free(w);
     958     
     959      angle = Projection(Tubevector[axis[1]], Vector[axis[0]]);
     960      fprintf(stdout, "Projection Tubevector1 axis[0] %lg %lg\n", angle, 1./angle);
     961      angle = Projection(Tubevector[axis[1]], Vector[axis[1]]);
     962      fprintf(stdout, "Projection Tubevector1 axis[1] %lg %lg\n", angle, 1./angle);
    887963     
    888964/*      fprintf(stdout, "Vector\n");
    … …  
    9881064      Tubevector[axis[2]][i] = Vector[axis[2]][i];
    9891065    }
    990    
     1066    // here we assume, that Vector[axis[2]] is along z direction!
     1067    gsl_matrix *M = gsl_matrix_alloc(2,2);
     1068    gsl_matrix *C = gsl_matrix_alloc(2,2);
     1069    gsl_matrix *evec = gsl_matrix_alloc(2,2);
     1070    gsl_vector *eval = gsl_vector_alloc(2);
     1071    gsl_vector *v = gsl_vector_alloc(2);
     1072    gsl_vector *u = gsl_vector_alloc(2);
     1073    gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(2);
     1074    gsl_matrix_set(C, 0,0, Vector[axis[0]][0]);
     1075    gsl_matrix_set(C, 1,0, Vector[axis[0]][1]);
     1076    gsl_matrix_set(C, 0,1, Vector[axis[1]][0]);
     1077    gsl_matrix_set(C, 1,1, Vector[axis[1]][1]);
     1078    gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0, C, C, 0.0, M);
     1079    fprintf(stdout, "M: \t%lg\t%lg\n\t%lg\t%lg\n", gsl_matrix_get(M,0,0), gsl_matrix_get(M,0,1), gsl_matrix_get(M,1,0), gsl_matrix_get(M,1,1));
     1080    gsl_eigen_symmv(M, eval, evec, w);
     1081    gsl_eigen_symmv_sort(eval,evec,GSL_EIGEN_SORT_ABS_DESC);
     1082    fprintf(stdout, "Eigenvalues: %lg\t%lg\n", gsl_vector_get(eval,0), gsl_vector_get(eval,1));
     1083    fprintf(stdout, "Eigenvectors: \t%lg\t%lg\n\t\t%lg\t%lg\n", gsl_matrix_get(evec,0,0), gsl_matrix_get(evec,0,1), gsl_matrix_get(evec,1,0), gsl_matrix_get(evec,1,1));
     1084    gsl_matrix_set(M, 0,0, 0.);
     1085    gsl_matrix_set(M, 1,0, 1.);
     1086    gsl_matrix_set(M, 0,1, -gsl_vector_get(eval,1)/gsl_vector_get(eval,0));
     1087    gsl_matrix_set(M, 1,1, 0.);
     1088    gsl_vector_set(v,0,(double)chiral[0]);
     1089    gsl_vector_set(v,1,(double)chiral[1]);
     1090    gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0, evec, M, 0.0, C);
     1091    gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0, C, evec, 0.0, M);
     1092    fprintf(stdout, "M: \t%lg\t%lg\n\t%lg\t%lg\n", gsl_matrix_get(M,0,0), gsl_matrix_get(M,0,1), gsl_matrix_get(M,1,0), gsl_matrix_get(M,1,1));
     1093    gsl_blas_dgemv(CblasNoTrans, 1.0, M, v, 0.0, u);
     1094    fprintf(stdout, "Looking for factor to integer...\n");
     1095    for(i=1;i<(chiral[0]+chiral[1])*(chiral[0]+chiral[1]);i++) {
     1096      x1 = gsl_vector_get(u,0)*(double)i;
     1097      x2 = gsl_vector_get(u,1)*(double)i;
     1098      x3 =
     1099      fprintf(stdout, "%d: %d\t%d vs. %lg\t%lg\n",i, ((int)(x1+x1/fabs(x1)*.5)), ((int)(x2+x2/fabs(x2)*.5)), (x1), (x2));
     1100      if (( fabs( ((int)(x1+x1/fabs(x1)*.5)) - (x1) ) < 1e-3) && ( fabs( ((int)(x2+x2/fabs(x2)*.5)) - (x2) ) < 1e-3 )) {
     1101        gsl_blas_dscal((double)i, u);
     1102        break;
     1103      }
     1104    }
     1105    fprintf(stdout, "(c,d) = (%lg,%lg)\n",gsl_vector_get(u,0), gsl_vector_get(u,1));
     1106
     1107    // get length
     1108    double x[NDIM];
     1109    for (i=0;i<NDIM;i++)
     1110      x[i] = gsl_vector_get(u,0) * Vector[axis[0]][i] + gsl_vector_get(u,1) * Vector[axis[1]][i];
     1111    angle = Norm(x)/Norm(Tubevector[axis[1]]) ;//ScalarProduct(x,Tubevector[axis[1]])/Norm(Tubevector[axis[1]]);
     1112    fprintf(stdout, "angle is %lg\n", angle);
     1113    for (i=0;i<NDIM;i++) {
     1114      Tubevector[axis[1]][i] = gsl_vector_get(u,0) * Vector[axis[0]][i] + gsl_vector_get(u,1) * Vector[axis[1]][i];
     1115    }
     1116
     1117    // Probe
     1118    gsl_matrix_set(M, 0,0, Vector[axis[0]][0]);
     1119    gsl_matrix_set(M, 1,0, Vector[axis[0]][1]);
     1120    gsl_matrix_set(M, 0,1, Vector[axis[1]][0]);
     1121    gsl_matrix_set(M, 1,1, Vector[axis[1]][1]);
     1122    gsl_vector_set(v,0,(double)chiral[0]);
     1123    gsl_vector_set(v,1,(double)chiral[1]);
     1124    gsl_blas_dgemv(CblasNoTrans, 1.0, M, u, 0.0, eval);
     1125    gsl_blas_dgemv(CblasTrans, 1.0, M, v, 0.0, u);
     1126    x1=1.;
     1127    gsl_blas_ddot(u,eval,&x1);
     1128    fprintf(stdout, "Testing (c,d): (a,b) M^t M (c,d)^t = 0 ? : %lg\n", x1);
     1129
     1130    gsl_matrix_free(M);
     1131    gsl_matrix_free(C);
     1132    gsl_matrix_free(evec);
     1133    gsl_vector_free(eval);
     1134    gsl_vector_free(v);
     1135    gsl_vector_free(u);
     1136    gsl_eigen_symmv_free(w);   
    9911137
    9921138    // retrieve seed ...
    … …  
    11511297            //PrintVector(stdout, coord);
    11521298            // project down on major and minor Tubevectors and check for length if out of sheet
    1153             tempvector = MatrixTrafoInverse(offset, TubevectorInverse);
    1154             if (((tempvector[axis[0]] + MYEPSILON) >= 0) && ((factors[0] - tempvector[axis[0]]) > MYEPSILON) &&
    1155                 ((tempvector[axis[1]] + MYEPSILON) >= 0) && ((factors[1] - tempvector[axis[1]]) > MYEPSILON) &&
    1156                 ((tempvector[axis[2]] + MYEPSILON) >= 0) && ((factors[2] - tempvector[axis[2]]) > MYEPSILON)) { // check if within rotated cell          numbersheet++;
     1299            tempvector = MatrixTrafoInverse(coord, TubevectorInverse);
     1300            if (((tempvector[axis[0]] + MYEPSILON) > 0) && ((factors[0] - tempvector[axis[0]]) > MYEPSILON) &&
     1301                ((tempvector[axis[1]] + MYEPSILON) > 0) && ((factors[1] - tempvector[axis[1]]) > MYEPSILON) &&
     1302                ((tempvector[axis[2]] + MYEPSILON) > 0) && ((factors[2] - tempvector[axis[2]]) > MYEPSILON)) { // check if within rotated cell          numbersheet++;
    11571303              //if (nummer >= 2)  strcpy(atombuffer[nr].name, "O");
    11581304              //nummer++;
Note: See TracChangeset for help on using the changeset viewer.