saturn ring

acn_funk_head.h

@@ -0,0 +1,93 @@
+//headers of all used functions
+
+//num_rec.c 
+int tqli(float d[], float e[], int n, float **z);
+void tred2(float **a, int n, float d[], float e[]);
+float pythag(float a, float b);
+float ran2(long *idum);
+
+//acn_funk_lb
+void zacPriblUF(long seed);
+void calcF2Ul(int l, float **f);
+void calcFNEW2F();
+void fillHl(int l,float (*h)[4]);
+void izrSStressL(int l, float (*sigl)[4], float (*qlm)[4], float *ul, float *dx, float *dy, float *dz);
+void izrAStressL(int l,float (*tau)[4]);
+void izrOdvAStressL(int l,float *divTaul);
+void izrFeq(float **feq);
+void izrP();
+void fillUl2(int l, float* Ul2);
+int calcLBlnew(int l, int m);
+int calcLBmnew(int m);
+int calcLBmnew_noslipZ(int m);
+int calcm2mrot(int m, int faktor);
+int calcLpbc(int l);
+void pbcLB();
+void izrCasKorU();
+void izrQU();
+void compute_stress_tensor();
+
+//acn_funk_op.c
+void zacPriblQ(long seed);
+void dir2QRP(int l, float dirX, float dirY, float dirZ, float nX, float nY, float nZ, float vX, float vY, float vZ, char lmark, long seed);
+void dir2QBulk(int l,long seed);
+void izrOdvQ2(int l, float* dxx, float* dyy, float* dzz, float* dx, float* dy, float* dz);
+// void izrOdvQ2_fe(int l, float* dxx, float* dyy, float* dzz, float* dxy, float* dxz, float* dyz, float* dx, float* dy, float* dz);
+void izrOdvQ2POV(int l, float* dx, float* dy, float* dz);
+double calculate_free_energy();
+void izrOdvQ2_lb(int l, float* dxx, float* dyy, float* dzz, float* dx, float* dy, float* dz);
+void izrOdvU2(int l, float *ux, float *uy, float *uz);
+void fillQl2(int l, float* ql);
+void fillQlm2(int l, float (*qlm)[4]);
+void calcU1(int l,float *u1, float* dxx, float* dyy, float* dzz, float* dx, float* dy, float* dz, float* ql);
+void calcU2(float *u2, float *ux, float *uy, float *uz, float (*qlm)[4]);
+void calcU3(float *u3, float* dx, float* dy, float* dz, float* ul);
+void izrCasKorQ();
+void calcQNEW2Q();
+
+//acn_utils.c
+int i_vr(int l);
+int j_vr(int l);
+int k_vr(int l);
+void initialiseE();
+void AMatrixB(float (*c)[4], float (*a)[4],float (*b)[4]);
+
+//acn_funk_zapis.c
+void zapisQ(char *ime);
+void zapisTEN(char *ime,float **q, int lmax, int mmax);
+void beriDATA(char *ime);
+void zapisDATA(char *ime);
+void za_profilU(char *ime, int i, int j);
+void zapisUraw(char *ime);
+void zapisDENraw(char *ime);
+void calcQ2dirL(int l,float *dir);
+void zapisQ2DIRraw(char *ime);
+void zapisQ2OPraw(char *ime);
+void zapisPOLJEraw(int t,char *ime);
+//acn_funk_pe.c
+
+
+//acn_funk_kol.c
+
+//active droplet
+struct active_droplet{
+    
+    int num_points;  // number of lattice points on sphere
+    int* association;
+    double* lattice_thetas;
+    double* lattice_phis;
+    double* velocities;
+    double defect_config[13];
+
+};
+
+typedef struct active_droplet ActiveDroplet;
+
+
+
+
+// //active droplet
+// void propagate_trajectories_block(double* variables, double nu, double chi, double xi_rot, double xi_alpha,
+//   double dt, int num_steps);
+// void velocity_field(double* defect_config, int n, double* lattice_thetas, double* lattice_phis,
+//   double rho, double* velocities);

diff_eq_system.h

@@ -0,0 +1,4 @@
+void propagate_trajectories_block(double* variables, double nu, double chi, double xi_rot, double xi_alpha,
+  double dt, int num_steps);
+void velocity_field(double* defect_config, int n, double* lattice_thetas, double* lattice_phis,
+  double rho, double* velocities);

dipol/fibonacci/saturn_ring/state_A1_test/acn_active_droplet.c

@@ -0,0 +1,608 @@
+#include "diff_eq_system.h"
+
+void setDefectConfig(double* config, double* defect_config) {
+    for (int i=0; i<13; i++) {
+    defect_config[i] = config[i];
+    }
+}
+
+
+//kartezicne koo.v sfericne
+void setLattice(int l, double* spherical_lattice, double* lattice_thetas, double* lattice_phis) {
+    double xy;
+    for (int i=0; i<l; i++) {
+    xy = sqrt(spherical_lattice[3*i] * spherical_lattice[3*i]
+                + spherical_lattice[3*i+1] * spherical_lattice[3*i+1]);
+    lattice_thetas[i] = atan2(xy, spherical_lattice[3*i+2]);
+    lattice_phis[i] = atan2(spherical_lattice[3*i+1], spherical_lattice[3*i]);
+    }
+}
+
+
+int count_sphere_points(double* center, double radius){
+    double distance;
+    int count_points = 0;
+//     #pragma omp parallel for reduction (+:count_points)
+    for(int l=0;l<NMAX;l++){
+        
+        float x,y,z,r; 
+            
+        x=(float)(i_vr(l))-(float)(I-1)/2+0.001;
+        y=(float)(j_vr(l))-(float)(J-1)/2+0.001;
+        z=(float)(k_vr(l))-(float)(K-1)/2+0.001;
+        
+        distance = (x - center[0]) * (x - center[0]) + (y - center[1]) * (y - center[1]) + (z - center[2]) * (z - center[2]);
+        
+        if(LMARK[l]&(LMARKPOV) && distance < 1.01 * radius * radius){count_points++;}
+        
+    }    
+    return count_points;
+}
+
+
+void points_sphere(double* center, double radius, double* points, int* association){
+
+    int k = 0;
+    for(int l=0;l<NMAX;l++){
+        
+        if(LMARK[l]&(LMARKPOV)){
+            
+            float x,y,z,r; 
+            
+            x=(float)(i_vr(l))-(float)(I-1)/2+0.001;
+            y=(float)(j_vr(l))-(float)(J-1)/2+0.001;
+            z=(float)(k_vr(l))-(float)(K-1)/2+0.001;
+            
+            if ((x - center[0]) * (x - center[0]) + (y - center[1]) * (y - center[1]) + (z - center[2]) * (z - center[2]) < 1.01 * radius * radius) {
+//                 LMARK[l]=LMARKSPH;
+                
+                points[3*k]=x;
+                points[3*k+1]=y;
+                points[3*k+2]=z;
+                association[k] = l;
+                k++;
+                
+                
+            }
+        
+    }
+    }
+    
+}
+
+
+
+void active_droplet_var(ActiveDropletVariables* active_variables){
+//     active_variables->
+    active_variables->active_nu = 2;
+    active_variables->active_tau = 2;
+    active_variables->drag_mu = 2;
+    active_variables->defect_size = 0.001;
+    active_variables->xi_rot = 4 * active_variables->active_tau / log(1 / active_variables->defect_size);
+    active_variables->xi_alpha = active_variables->xi_rot;
+    active_variables->active_chi = 2 / log(1 / active_variables->defect_size);    
+    active_variables->active_dt = dt_active0 * dt_active0 * active_variables->active_nu * R / VELOCITY_SCALING;  // timestep for active simulation for defect velocity to mach rotation time
+    
+}
+
+// everything on a sphere
+void init_droplet(ActiveDroplet* sphere) {
+    
+    // double center1[3] = {0., 0., 0.};
+    double original_config[13] = {2.2549091,
+                                    2.5599512,
+                                    0.8845974,
+                                    0.5845573,
+                                    0.9355829,
+                                    3.9472116,
+                                    4.7124571,
+                                    1.7150761,
+                                    0.8025151,
+                                    0.6296737,
+                                    -2.3309994,
+                                    -2.5228739,
+                                    2.6938891};
+    
+    
+    int num_points = count_sphere_points(sphere->center, sphere->radius);
+    double* all_points = (double *)malloc(3 * num_points*sizeof(double));
+    int* assoc = (int *)malloc(num_points*sizeof(int));
+    points_sphere(sphere->center, sphere->radius, all_points, assoc);
+
+    sphere->num_points = num_points;
+    sphere->association = (int*)malloc(num_points*sizeof(int));
+    for (int i=0;i<num_points; i++) {
+        sphere->association[i]=assoc[i];
+//         if (i<100) {
+//             printf("%d\n", sphere->association[i]);
+//         }
+    }
+    sphere->lattice_thetas = (double *)malloc(num_points*sizeof(double));
+    sphere->lattice_phis = (double *)malloc(num_points*sizeof(double));
+    sphere->velocities = (double *)malloc(3*num_points*sizeof(double));
+                                
+    setLattice(sphere->num_points, all_points, sphere->lattice_thetas, sphere->lattice_phis);
+    setDefectConfig(original_config, sphere->defect_config);
+    
+    sphere->lattice_basis = (double *)malloc(num_points*9*sizeof(double));
+    lattice_base(sphere->num_points, sphere->lattice_thetas, sphere->lattice_phis, sphere->lattice_basis);
+    
+    sphere->theta_grad = (double *)malloc(sphere->num_points * sizeof(double));
+    sphere->phi_grad = (double *)malloc(sphere->num_points * sizeof(double));
+
+}
+
+void velocity_grad(ActiveDroplet* sphere){
+
+    for(int l=0;l<NMAX;l++){
+        
+        gradU[1][l]=0;
+        gradU[2][l]=0;
+        gradU[3][l]=0;
+    }
+     
+    double grad_x;
+    double grad_y;
+    double grad_z;
+    double vec_theta[3];
+    double vec_phi[3];
+    
+    #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500)
+    for(int k=0;k<sphere->num_points;k++){
+        float ux[4],uy[4],uz[4];
+        float n1,n2,n3;
+        
+
+        float pos[4];
+        pos[1]=(float)(i_vr(sphere->association[k])); //x
+        pos[2]=(float)(j_vr(sphere->association[k])); //y
+        pos[3]=(float)(k_vr(sphere->association[k])); //z
+        float n[4];
+        n[1]=QNIC[4][sphere->association[k]];
+        n[2]=QNIC[5][sphere->association[k]];
+        n[3]=QNIC[6][sphere->association[k]];
+        
+        float grad[4];
+        
+        
+        izrOdvU2smeri(pos, n, grad);
+        
+
+        gradU[1][sphere->association[k]]=grad[1];
+        gradU[2][sphere->association[k]]=grad[2];
+        gradU[3][sphere->association[k]]=grad[3];
+        
+        etheta(sphere->lattice_phis[k], sphere->lattice_thetas[k], vec_theta);
+        ephi(sphere->lattice_phis[k], vec_phi);
+        
+        
+        // Calculated gradient has direction opposite to flow velocity. In propagation code it should be taken in the direction of flow, so we add a minus.
+        // Gradients are also scaled with the velocity_ratio factor and DT to transform magnitudes to appropriate units.
+        sphere->theta_grad[k] = -(grad[1]*vec_theta[0] + grad[2]*vec_theta[1] + grad[3]*vec_theta[2])/VELOCITY_SCALING*DT;
+        sphere->phi_grad[k] = -(grad[1]*vec_phi[0] + grad[2]*vec_phi[1] + grad[3]*vec_phi[2])/VELOCITY_SCALING*DT;      
+        
+//         fprintf(stderr, "gradient %f\n",grad[1]);
+//         fprintf(stderr, "lokacija %d\n",(int)(pos[1]));
+//         fprintf(stderr,"%f\n", (grad[1]*n[1]+grad[2]*n[2]+grad[3]*n[3])/sqrt(grad[1]*grad[1]+grad[2]*grad[2]+grad[3]*grad[3]));
+        
+        
+
+    };
+
+}
+
+void defect_grad(ActiveDroplet* sphere){
+
+    double vec_theta[3];
+    double vec_phi[3];
+
+    float pos_defect[4];
+    FILE *defect_con;
+    FILE *defect_all_grad;
+
+    for(int j=0;j<4;j++){
+        double loc_defect[3];
+        double r=sphere->radius;
+        
+        sph_to_cart(sphere->defect_config[j], sphere->defect_config[j+4], loc_defect);            
+
+        pos_defect[1]=loc_defect[0]*r + (float)(I-1)/2+0.001; //x
+        pos_defect[2]=loc_defect[1]*r + (float)(J-1)/2+0.001; //y
+        pos_defect[3]=loc_defect[2]*r + (float)(K-1)/2+0.001; //z
+   
+        float n[4];
+        n[1]=loc_defect[0];
+        n[2]=loc_defect[1];
+        n[3]=loc_defect[2];
+        
+        float grad_defect[4];
+        izrOdvU2smeri(pos_defect, n, grad_defect);
+        
+        
+        etheta(sphere->defect_config[j+4], sphere->defect_config[j], vec_theta);
+        ephi(sphere->defect_config[j+4], vec_phi);
+        
+        // Calculated gradient has direction opposite to flow velocity. In velocity field code it should be taken in the direction of flow, so we add a minus.
+        // Gradients are also scaled with the velocity_ratio factor and DT to transform magnitudes to appropriate units.
+        sphere->defect_theta_grad[j] = -(grad_defect[1]*vec_theta[0] + grad_defect[2]*vec_theta[1] + grad_defect[3]*vec_theta[2])/VELOCITY_SCALING*DT;
+        sphere->defect_phi_grad[j] = -(grad_defect[1]*vec_phi[0] + grad_defect[2]*vec_phi[1] + grad_defect[3]*vec_phi[2])/VELOCITY_SCALING*DT;     
+        
+        sphere->defect_gradient[j] = sqrt(grad_defect[1]*grad_defect[1]+grad_defect[2]*grad_defect[2]+grad_defect[3]*grad_defect[3])/VELOCITY_SCALING*DT;
+
+    }
+
+}
+
+void zapis_all_grad(ActiveDroplet* sphere){
+
+    float pos_defect[4];
+    FILE *defect_con;
+    FILE *defect_all_grad;
+
+    float all_grad_defect[12];
+    
+    for(int j=0;j<4;j++){
+        double loc_defect[3];
+        double r=sphere->radius;
+        
+        sph_to_cart(sphere->defect_config[j], sphere->defect_config[j+4], loc_defect);            
+        
+        pos_defect[1]=loc_defect[0]*r + (float)(I-1)/2+0.001; //x
+        pos_defect[2]=loc_defect[1]*r + (float)(J-1)/2+0.001; //y
+        pos_defect[3]=loc_defect[2]*r + (float)(K-1)/2+0.001; //z    
+        
+        float n[4];
+        n[1]=loc_defect[0];
+        n[2]=loc_defect[1];
+        n[3]=loc_defect[2];
+        
+        float grad_defect[4];
+        
+        
+        izrOdvU2smeri(pos_defect, n, grad_defect);
+        
+        all_grad_defect[3*j] = grad_defect[1]/VELOCITY_SCALING*DT;
+        all_grad_defect[3*j+1] = grad_defect[2]/VELOCITY_SCALING*DT;
+        all_grad_defect[3*j+2] = grad_defect[3]/VELOCITY_SCALING*DT;
+
+        
+    }
+    
+    #if TWRITEgrad
+      //         zapis defect_config;
+        defect_con=fopen("./data/defect_config.txt","a");
+        
+        if(defect_con==NULL){
+            fprintf(stderr, "Failed opening file. a mode. file=%s\n", "./defect_con.txt");
+            exit(1); 
+        }
+        for(int c=0;c<13;c++){
+            fprintf(defect_con, "%lf;",sphere->defect_config[c]);
+        };
+        
+        fprintf(defect_con, "\n");
+        fclose(defect_con);
+        
+      //         zapis defect_all_grad;  
+        defect_all_grad=fopen("./data/defect_all_grad.txt","a");
+        if(defect_all_grad==NULL){
+            fprintf(stderr, "Failed opening file. a mode. file=%s\n", "./defect_all_grad.txt");
+            exit(1); 
+        }
+        for(int c=0;c<12;c++){
+            fprintf(defect_all_grad, "%lf;",all_grad_defect[c]);
+        };
+        
+//         fprintf(defect_all_grad, "%f;",all_grad_defect);
+        fprintf(defect_all_grad, "\n");
+        fclose(defect_all_grad);
+
+       
+    #endif
+
+}
+
+void velocity_sphere(ActiveDroplet* sphere, ActiveDropletVariables* active_variables){
+//     propagate_trajectories_block(sphere->defect_config, nu, chi, xi_rot, xi_alpha, dt, num_steps);
+//     double active_nu = 1.5;
+//     double active_tau = 2;
+//     double drag_mu = 2;
+//     double defect_size = 0.001;
+//     double xi_rot = 4 * active_tau / log(1 / defect_size);
+//     double xi_alpha = xi_rot;
+//     double active_chi = 2 / log(1 / defect_size);
+//     
+//     double active_dt = dt_active0*dt_active0*active_nu*R/VELOCITY_SCALING;  // timestep for active simulation for defect velocity to mach rotation time
+//     
+//     printf("REAL xi_rot = %f\n", xi_rot);
+ 
+    
+    double gradient=0;
+    for(int j=0;j<4;j++){gradient+=sphere->defect_gradient[j];}
+    
+    double active_nu_eff = active_variables->active_nu - active_variables->drag_mu * gradient/4;
+    
+//     fprintf(stderr, "active_nu_eff %f\n",active_nu_eff);
+    defect_grad(sphere);
+    
+    propagate_trajectories_step_drag(sphere->defect_config, sphere->defect_theta_grad, sphere->defect_phi_grad,
+                                     active_variables->active_nu, active_variables->active_chi, active_variables->drag_mu, 
+                                     active_variables->xi_rot, active_variables->xi_alpha, active_variables->active_dt);
+ 
+        
+    
+    // propagate_trajectories_block(sphere->defect_config, active_nu, active_chi, xi_rot, xi_alpha, dt, num_steps);
+    
+    velocity_grad(sphere);
+    
+    velocity_field_drag(sphere->defect_config, sphere->num_points, sphere->lattice_thetas, sphere->lattice_phis, 
+                        sphere->lattice_basis, sphere->theta_grad, sphere->phi_grad, 
+                        active_variables->active_nu, active_variables->drag_mu, 0.1, sphere->velocities, STPROC); //rho=0.1
+    
+//         velocity_field_fast(sphere->defect_config, sphere->num_points, sphere->lattice_thetas, sphere->lattice_phis,
+//                             sphere->lattice_basis, active_nu, 0.1, sphere->velocities, STPROC); //rho=0.1
+                            
+    // velocity_field(sphere->defect_config, sphere->num_points, sphere->lattice_thetas, sphere->lattice_phis,active_nu, 0.1, sphere->velocities); 
+
+    
+//     x_to_z_transf(sphere->num_points, sphere->velocities);
+    
+//     #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500)
+    for(int k=0;k<sphere->num_points;k++){   
+            float n1=QNIC[4][sphere->association[k]];	//Surface normal
+            float n2=QNIC[5][sphere->association[k]];
+            float n3=QNIC[6][sphere->association[k]];
+
+            UNIC[1][sphere->association[k]]=VELOCITY_SCALING*sphere->velocities[3*k];		//surface velocity
+            UNIC[2][sphere->association[k]]=VELOCITY_SCALING*sphere->velocities[3*k+1];
+            UNIC[3][sphere->association[k]]=VELOCITY_SCALING*sphere->velocities[3*k+2];
+            
+            
+//             if (k<1000) {printf("%f\n", sphere->velocities[3*k]*n1+sphere->velocities[3*k+1]*n2+sphere->velocities[3*k+2]*n3);}
+
+//             if (k==1000) {printf("%f ", sphere->velocities[3*k]*n1+sphere->velocities[3*k+1]*n2+sphere->velocities[3*k+2]*n3);}
+            
+
+
+    }
+    
+    
+}
+
+
+//Writes the velocity GRADIENT (with Paraview header)
+void zapisUgrad(char *ime){
+  // do it in binary
+  int l;
+  float grad[3];
+  FILE *file;
+
+  file=fopen(ime,"w");
+  if(file==NULL){
+    fprintf(stderr, "Failed opening file. w mode. file=%s\n", ime);
+     exit(1);
+  }
+
+  
+  fprintf(file,"# vtk DataFile Version 3.0\nKnot\nBINARY\nDATASET STRUCTURED_POINTS\n");
+  fprintf(file,"DIMENSIONS %d %d %d\n",I,J,K);
+  fprintf(file,"ORIGIN 0 0 0\n");
+  fprintf(file,"SPACING 1 1 1\n");
+  fprintf(file,"POINT_DATA %d\n",I*J*K);
+  fprintf(file,"VECTORS gradU float\n");
+  
+  
+  for(l=0; l<NMAX; l++){
+    // go to bigendian
+    grad[0] = FloatSwap(gradU[1][l]); 
+    grad[1] = FloatSwap(gradU[2][l]);
+    grad[2] = FloatSwap(gradU[3][l]);
+      
+    
+    fwrite(grad,sizeof(float),3,file);
+  };
+  fclose(file);
+}
+
+//Writes the velocity GRADIENT - DEFECT (with Paraview header)
+void zapisUgradDEFECT(char *ime){
+  // do it in binary
+  int l;
+  float grad_defect[3];
+  FILE *file;
+
+  file=fopen(ime,"w");
+  if(file==NULL){
+    fprintf(stderr, "Failed opening file. w mode. file=%s\n", ime);
+     exit(1);
+  }
+
+  
+  fprintf(file,"# vtk DataFile Version 3.0\nKnot\nBINARY\nDATASET STRUCTURED_POINTS\n");
+  fprintf(file,"DIMENSIONS %d %d %d\n",I,J,K);
+  fprintf(file,"ORIGIN 0 0 0\n");
+  fprintf(file,"SPACING 1 1 1\n");
+  fprintf(file,"POINT_DATA %d\n",I*J*K);
+  fprintf(file,"VECTORS gradUdefect float\n");
+  
+
+  for(l=0; l<NMAX; l++){
+    // go to bigendian
+    grad_defect[0] = FloatSwap(gradU_defect[1][l]); 
+    grad_defect[1] = FloatSwap(gradU_defect[2][l]);
+    grad_defect[2] = FloatSwap(gradU_defect[3][l]);
+      
+    
+    fwrite(grad_defect,sizeof(float),3,file);
+  };
+  fclose(file);
+}
+
+
+
+void locDefect(ActiveDroplet* sphere, char *ime){
+    FILE *file;
+    
+    file=fopen(ime,"w");
+    if(file==NULL){
+        fprintf(stderr, "Failed opening file. w mode. file=%s\n", ime);
+        exit(1);
+    }
+    
+    fprintf(file,"x,y,z\n");
+    //4 defects
+    for(int i=0;i<4;i++){
+        double loc_defect[3];
+        double r=sphere->radius;
+        
+        sph_to_cart(sphere->defect_config[i], sphere->defect_config[i+4], loc_defect);    
+        fprintf(file,"%f,%f,%f\n",loc_defect[0]*r + sphere->center[0], loc_defect[1]*r + sphere->center[1], loc_defect[2]*r + sphere->center[2]);
+    }
+    fclose(file);
+    
+    
+}
+
+
+//----------------------------------------FIBONACCI------------------------------------------------------------------------
+
+void fibonacci_init(FibonacciLattice* fibonacci_lattice, int num_points, double center[3], double radius){
+    double golden_ratio = 1.618033988749894;
+    fibonacci_lattice->center = center;
+    fibonacci_lattice->radius = radius;
+    fibonacci_lattice->num_points = num_points;
+    
+    fibonacci_lattice->lattice_thetas = (double *)malloc(num_points*sizeof(double));
+    fibonacci_lattice->lattice_phis = (double *)malloc(num_points*sizeof(double));
+    
+    for(int i=0;i<num_points;i++){
+        fibonacci_lattice->lattice_phis[i] = (i+1)*(2-golden_ratio)*2*M_PI;
+        fibonacci_lattice->lattice_thetas[i] = acos(2*(i+1)/(num_points+1)-1);
+        
+    }
+    
+    fibonacci_lattice->velocities = (double *)malloc(3*num_points*sizeof(double));
+    
+    fibonacci_lattice->lattice_basis = (double *)malloc(num_points*9*sizeof(double));
+    lattice_base(fibonacci_lattice->num_points, fibonacci_lattice->lattice_thetas, fibonacci_lattice->lattice_phis, fibonacci_lattice->lattice_basis);
+    
+    fibonacci_lattice->theta_grad = (double *)malloc(fibonacci_lattice->num_points * sizeof(double));
+    fibonacci_lattice->phi_grad = (double *)malloc(fibonacci_lattice->num_points * sizeof(double));
+}
+
+
+void grad_fibonacci(FibonacciLattice* fibonacci_lattice){
+    double vec_theta[3];
+    double vec_phi[3];
+    float position[4];
+    FILE *defect_con;
+    FILE *defect_all_grad;
+/*    
+    for(int l=0;l<NMAX;l++){
+        
+        gradU_defect[1][l]=0;
+        gradU_defect[2][l]=0;
+        gradU_defect[3][l]=0;
+    }*/
+    
+    for(int j=0;j<fibonacci_lattice->num_points;j++){
+        
+        double location[3]; //normirane kartezicne koordinate
+        double r=fibonacci_lattice->radius;
+        
+        sph_to_cart(fibonacci_lattice->lattice_thetas[j], fibonacci_lattice->lattice_phis[j], location);            
+        
+        position[1]=location[0]*r + (float)(I-1)/2+0.001; //x
+        position[2]=location[1]*r + (float)(J-1)/2+0.001; //y
+        position[3]=location[2]*r + (float)(K-1)/2+0.001; //z
+/*        fprintf(stderr, "position = %f, %f, %f;   ",position[1],position[2],position[3]);
+        fprintf(stderr, "%f\n ", position[1]*position[1]+position[2]*position[2]+position[3]*position[3]);
+ */      
+        float n[4];
+        n[1]=location[0];
+        n[2]=location[1];
+        n[3]=location[2];
+        
+        float grad_point[4];
+        izrOdvU2smeri(position, n, grad_point);
+
+        
+        for(int i=0;i<3;i++){
+            vec_theta[i] = fibonacci_lattice->lattice_basis[9*j+3+i];
+            vec_phi[i] = fibonacci_lattice->lattice_basis[9*j+6+i];
+        }
+        
+
+        
+        // Calculated gradient has direction opposite to flow velocity. In velocity field code it should be taken in the direction of flow, so we add a minus.
+        // Gradients are also scaled with the velocity_ratio factor and DT to transform magnitudes to appropriate units.
+        fibonacci_lattice->theta_grad[j] = -(grad_point[1]*vec_theta[0] + grad_point[2]*vec_theta[1] + grad_point[3]*vec_theta[2])/VELOCITY_SCALING*DT;
+        fibonacci_lattice->phi_grad[j] = -(grad_point[1]*vec_phi[0] + grad_point[2]*vec_phi[1] + grad_point[3]*vec_phi[2])/VELOCITY_SCALING*DT;     
+           
+    }
+
+}
+
+
+
+void velocity_fibonacci(ActiveDroplet* sphere, FibonacciLattice* fibonacci_lattice, ActiveDropletVariables* active_variables){
+    
+    grad_fibonacci(fibonacci_lattice);
+
+    
+    velocity_field_drag(sphere->defect_config, fibonacci_lattice->num_points, fibonacci_lattice->lattice_thetas, fibonacci_lattice->lattice_phis, 
+                        fibonacci_lattice->lattice_basis, fibonacci_lattice->theta_grad, fibonacci_lattice->phi_grad, 
+                        active_variables->active_nu, active_variables->drag_mu, 0.1, fibonacci_lattice->velocities, STPROC); //rho=0.1
+}
+
+
+//average fibonacci_lattice->velocities //vector 
+void zapis_fibonacci_average_velocity(ActiveDroplet* sphere, FibonacciLattice* fibonacci_lattice, ActiveDropletVariables* active_variables){
+    FILE *file;
+    
+    float Vx = 0;
+    float Vy = 0;
+    float Vz = 0;
+    
+//     fprintf(stderr,"active_nu = %f\n", active_variables->active_nu);
+//     fprintf(stderr,"drag_mu = %f\n", active_variables->drag_mu);
+//     fprintf(stderr,"active_tau = %f\n", active_variables->active_tau);
+//     fprintf(stderr,"defect_size = %f\n", active_variables->defect_size);
+//     fprintf(stderr,"xi_rot = %f\n", active_variables->xi_rot);
+//     fprintf(stderr,"xi_alpha = %f\n", active_variables->xi_alpha);
+//     fprintf(stderr,"active_chi = %f\n", active_variables->active_chi);
+//     fprintf(stderr,"active_dt = %f\n", active_variables->active_dt);   
+//  
+   
+    velocity_fibonacci(sphere, fibonacci_lattice, active_variables);
+    
+    for(int i=0;i<fibonacci_lattice->num_points;i++){
+//         fprintf(stderr,"%f; ", fibonacci_lattice->velocities[3*i]);
+//         fprintf(stderr,"Vx = %f\n", Vx);
+        Vx = Vx + fibonacci_lattice->velocities[3*i];
+        Vy = Vy + fibonacci_lattice->velocities[3*i+1];
+        Vz = Vz + fibonacci_lattice->velocities[3*i+2];
+        
+    }
+    /*
+    fprintf(stderr,"%f,%f,%f\n", fibonacci_lattice->velocities[0],fibonacci_lattice->velocities[1],fibonacci_lattice->velocities[2]);
+        */
+/*    fprintf(stderr,"Vx = %f\n", Vx);
+    fprintf(stderr,"%f,%f,%f\n", Vx/fibonacci_lattice->num_points,Vy/fibonacci_lattice->num_points, Vz/fibonacci_lattice->num_points);
+ */       
+    
+    file=fopen("./data/fibonacci_velocity.txt","a");
+        
+    if(file==NULL){
+        fprintf(stderr, "Failed opening file. a mode. file=%s\n", "./fibonacci_velocity.txt");
+        exit(1); 
+    }
+    
+    fprintf(file, "%f,%f,%f\n",Vx/fibonacci_lattice->num_points, Vy/fibonacci_lattice->num_points, Vz/fibonacci_lattice->num_points);
+    
+    fclose(file);
+    
+    
+    
+    
+}

dipol/fibonacci/saturn_ring/state_A1_test/acn_funk_lb.c

@@ -0,0 +1,607 @@
+//Initialisation for the velocity field and partial distribution functions
+void zacPriblUF(long seed){
+  int l,m;
+  float perturbacija;
+  ran2(&seed);
+  for(l=0;l<NMAX;l++){   
+      float x,y,z, dirX, dirY, dirZ;
+  
+      x=(float)(i_vr(l))+(float)(I)/2;
+      y=(float)(j_vr(l))+(float)(J)/2;
+      z=(float)(k_vr(l))-4;
+  
+      dirX=0;//.01*x/sqrt(x*x+y*y+z*z);
+      dirY=0;//.01*y/sqrt(x*x+y*y+z*z);
+      dirZ=0;//.01*z/sqrt(x*x+y*y+z*z);
+    U[1][l]=0;//.1*ran2(&seed);//05f*( J*J/4-(j_vr(l)-J/2)*(j_vr(l)-J/2))/J/J;//ran2(&seed);
+    U[2][l]=0;//.1*ran2(&seed);//0.005f*ran2(&seed);
+    U[3][l]=0;//.1*ran2(&seed);
+    U[4][l]=DENSITYINIT;
+  }
+    
+   izrFeq(FEQ);
+   #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500) private(m)
+   for(l=0;l<NMAX;l++){
+    for(m=0;m<ST_HITROSTI;m++){F[m][l]=FEQ[m][l]; FNEW[m][l]=F[m][l]; }
+   }
+}
+
+
+//Computes velocity field from F
+void calcF2Ul(int l, float **f){
+  int m;
+  float fex=0.f, fey=0.f, fez=0.f, density=0.f;
+  
+  if(LMARK[l]==LMARK1){
+    for(m=0;m<ST_HITROSTI;m++){
+      density+=f[m][l];
+      fex+=f[m][l]*E[m][1];
+      fey+=f[m][l]*E[m][2];
+      fez+=f[m][l]*E[m][3];
+    }
+    float faktor=0.5f;	//Correction due to discrete lattice effects
+    fex+=faktor*FORCE[1][l];
+    fey+=faktor*FORCE[2][l];
+    fez+=faktor*FORCE[3][l];
+    
+    U[1][l]=fex/density;
+    U[2][l]=fey/density;
+    U[3][l]=fez/density;
+    U[4][l]=density;
+//     printf("%d ",LMARK[l]);
+
+  }
+ else{ U[4][l]=DENSITYINIT; U[1][l]=UNIC[1][l]; U[2][l]=UNIC[2][l]; U[3][l]=UNIC[3][l];}
+}
+
+
+void calcFNEW2F(){
+  int l;
+  #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500)
+  for(int l=0; l<NMAX; l++){
+    if(LMARK[l]==LMARK1) for(int m=0; m<ST_HITROSTI; m++) F[m][l]=FNEW[m][l];
+ }
+}
+
+
+//Fills in local molecular field from the global matrix
+void fillHl(int l,float (*h)[4])
+{
+	h[1][1]=H[1][l];
+	h[2][2]=H[2][l];
+	h[3][3]=H[3][l];
+	h[1][2]=H[4][l];
+	h[1][3]=H[5][l];
+	h[2][3]=H[6][l];
+	h[2][1]=h[1][2];
+	h[3][1]=h[1][3];
+	h[3][2]=h[2][3];
+}
+
+void compute_stress_tensor(){
+  #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500)
+  for(int l=0;l<NMAX;l++){
+    if(!(LMARK[l]==LMARK1)) continue;
+    float trQH;
+    float hl[4][4],qlm[4][4];
+    float qlm_hl[4][4];
+  //   float hl_qlm[4][4];
+
+    fillHl(l,hl);					//fill the local molecular field matrix
+    fillQlm2(l,qlm);				//fill the local order parameter tensor matrix
+
+    AMatrixB(qlm_hl,qlm,hl);		//Matrix multiplication
+  //   AMatrixB(hl_qlm,hl,qlm);
+    
+    trQH=qlm_hl[1][1]+qlm_hl[2][2]+qlm_hl[3][3];
+    
+    float tau_sim12=(-XI*(qlm_hl[2][1]+qlm_hl[1][2]+2.f/3.f*hl[1][2])+2.f*XI*(qlm[1][2])*trQH );
+    float tau_asim12=qlm_hl[1][2]-qlm_hl[2][1];
+    float tau_sim13=(-XI*(qlm_hl[3][1]+qlm_hl[1][3]+2.f/3.f*hl[1][3])+2.f*XI*(qlm[1][3])*trQH );
+    float tau_asim13=qlm_hl[1][3]-qlm_hl[3][1];
+    float tau_sim23=(-XI*(qlm_hl[3][2]+qlm_hl[2][3]+2.f/3.f*hl[2][3])+2.f*XI*(qlm[2][3])*trQH );
+    float tau_asim23=qlm_hl[2][3]-qlm_hl[3][2];
+
+
+    TAU[1][l]=(-XI*(2*qlm_hl[1][1]+2.f/3.f*hl[1][1])+2.f*XI*(qlm[1][1]+1/3.f)*trQH ) -ZETA*qlm[1][1];				//calculate local stress tensor
+    TAU[2][l]=(-XI*(2*qlm_hl[2][2]+2.f/3.f*hl[2][2])+2.f*XI*(qlm[2][2]+1/3.f)*trQH ) -ZETA*qlm[2][2];
+    TAU[3][l]=(-XI*(2*qlm_hl[3][3]+2.f/3.f*hl[3][3])+2.f*XI*(qlm[3][3]+1/3.f)*trQH ) -ZETA*qlm[3][3];
+    TAU[4][l]=tau_asim12+tau_sim12 -ZETA*qlm[1][2];
+    TAU[5][l]=tau_asim13+tau_sim13 -ZETA*qlm[1][3];
+    TAU[6][l]=tau_asim23+tau_sim23 -ZETA*qlm[2][3];
+    TAU[7][l]=-tau_asim12+tau_sim12 -ZETA*qlm[2][1];
+    TAU[8][l]=-tau_asim13+tau_sim13 -ZETA*qlm[3][1];
+    TAU[9][l]=-tau_asim23+tau_sim23 -ZETA*qlm[3][2];
+  }
+  
+  #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500)
+  for(int l=0; l<NMAX; l++){
+    int lpbc=calcLpbc(l);
+    if(lpbc!=l) for(int m=1; m<=9; m++) TAU[m][l]=TAU[m][lpbc];
+  }
+}
+
+//Funkcija, ki izracuna stress tenzor v dani tocki
+// void izrAStressL(int l,float (*tau)[4]){
+//   float trQH;
+//   float hl[4][4],qlm[4][4];
+//   float qlm_hl[4][4];
+//   float hl_qlm[4][4];
+// 
+//   fillHl(l,hl);					//fill the local molecular field matrix
+//   fillQlm2(l,qlm);				//fill the local order parameter tensor matrix
+// 
+//   AMatrixB(qlm_hl,qlm,hl);		//Matrix multiplication
+//   AMatrixB(hl_qlm,hl,qlm);
+//   
+//   trQH=qlm_hl[1][1]+qlm_hl[2][2]+qlm_hl[3][3];
+// 
+// 
+//   tau[1][1]=qlm_hl[1][1]-hl_qlm[1][1]+(-XI*(hl_qlm[1][1]+qlm_hl[1][1]+2.f/3.f*hl[1][1])+2.f*XI*(qlm[1][1]+1/3.f)*trQH ) -ZETA*qlm[1][1];				//calculate local stress tensor
+//   tau[2][2]=qlm_hl[2][2]-hl_qlm[2][2]+(-XI*(hl_qlm[2][2]+qlm_hl[2][2]+2.f/3.f*hl[2][2])+2.f*XI*(qlm[2][2]+1/3.f)*trQH ) -ZETA*qlm[2][2];
+//   tau[3][3]=qlm_hl[3][3]-hl_qlm[3][3]+(-XI*(hl_qlm[3][3]+qlm_hl[3][3]+2.f/3.f*hl[3][3])+2.f*XI*(qlm[3][3]+1/3.f)*trQH ) -ZETA*qlm[3][3];
+//   tau[1][2]=qlm_hl[1][2]-hl_qlm[1][2]+(-XI*(hl_qlm[1][2]+qlm_hl[1][2]+2.f/3.f*hl[1][2])+2.f*XI*(qlm[1][2])*trQH ) -ZETA*qlm[1][2];
+//   tau[1][3]=qlm_hl[1][3]-hl_qlm[1][3]+(-XI*(hl_qlm[1][3]+qlm_hl[1][3]+2.f/3.f*hl[1][3])+2.f*XI*(qlm[1][3])*trQH ) -ZETA*qlm[1][3];
+//   tau[2][3]=qlm_hl[2][3]-hl_qlm[2][3]+(-XI*(hl_qlm[2][3]+qlm_hl[2][3]+2.f/3.f*hl[2][3])+2.f*XI*(qlm[2][3])*trQH ) -ZETA*qlm[2][3];
+//   tau[2][1]=qlm_hl[2][1]-hl_qlm[2][1]+(-XI*(hl_qlm[1][2]+qlm_hl[1][2]+2.f/3.f*hl[1][2])+2.f*XI*(qlm[1][2])*trQH ) -ZETA*qlm[2][1];
+//   tau[3][1]=qlm_hl[3][1]-hl_qlm[3][1]+(-XI*(hl_qlm[1][3]+qlm_hl[1][3]+2.f/3.f*hl[1][3])+2.f*XI*(qlm[1][3])*trQH ) -ZETA*qlm[3][1];
+//   tau[3][2]=qlm_hl[3][2]-hl_qlm[3][2]+(-XI*(hl_qlm[2][3]+qlm_hl[2][3]+2.f/3.f*hl[2][3])+2.f*XI*(qlm[2][3])*trQH ) -ZETA*qlm[3][2];
+// }
+
+void izrOdvAStressL(int l,float *divTaul){
+  
+  
+  int lp,lm;
+  float fak;
+
+//   float taup[4][4],taum[4][4];
+  float DTAUl[4][4];
+  
+  fak=0.5f;	//TAU needs to be calculated only for LMARK1 markers
+
+  lp=l+1; lm=l-1;	
+  if(!(LMARK[lm]&(LMARK1))) {lm=l,fak=1.f;} else if(!(LMARK[lp]&LMARK1)) {lp=l;fak=1.f;};
+//   izrAStressL(lp,taup);									//calculate stress in lxp nad lxm points
+//   izrAStressL(lm,taum);
+  DTAUl[1][1]=fak*(TAU[1][lp]-TAU[1][lm]);							//calculate derivatives of stress point: eg DTAU[2][1]=DTau_yx /dx
+  DTAUl[2][1]=fak*(TAU[7][lp]-TAU[7][lm]);
+  DTAUl[3][1]=fak*(TAU[8][lp]-TAU[8][lm]);
+
+  lp=l+I; lm=l-I;
+  if(!(LMARK[lm]&LMARK1)) {lm=l,fak=1.f;} else if(!(LMARK[lp]&LMARK1)) {lp=l;fak=1.f;};
+//   izrAStressL(lp,taup);									//calculate stress in lyp nad lym points
+//   izrAStressL(lm,taum);
+  DTAUl[1][2]=fak*(TAU[4][lp]-TAU[4][lm]);							//calculate derivatives of stress point: eg DTAU[2][1]=DTau_yx /dx
+  DTAUl[2][2]=fak*(TAU[2][lp]-TAU[2][lm]);
+  DTAUl[3][2]=fak*(TAU[9][lp]-TAU[9][lm]);
+
+  lp=l+I*J; lm=l-I*J;
+  if(!(LMARK[lm]&LMARK1)) {lm=l,fak=1.f;} else if(!(LMARK[lp]&LMARK1)) {lp=l;fak=1.f;};
+//   izrAStressL(lp,taup);									//calculate stress in lyp nad lym points
+//   izrAStressL(lm,taum);
+  DTAUl[1][3]=fak*(TAU[5][lp]-TAU[5][lm]);							//calculate derivatives of stress point: eg DTAU[2][1]=DTau_yx /dx
+  DTAUl[2][3]=fak*(TAU[6][lp]-TAU[6][lm]);
+  DTAUl[3][3]=fak*(TAU[3][lp]-TAU[3][lm]);
+
+  divTaul[1]=DTAUl[1][1]+DTAUl[1][2]+DTAUl[1][3];
+  divTaul[2]=DTAUl[2][1]+DTAUl[2][2]+DTAUl[2][3];
+  divTaul[3]=DTAUl[3][1]+DTAUl[3][2]+DTAUl[3][3];
+}
+
+//Funkcija, ki izracuna odvode stress tenzorja
+// void izrOdvAStressL(int l,float *divTaul){
+//   
+//   
+//   int lp,lm;
+//   float fak;
+// 
+//   float taup[4][4],taum[4][4], DTAUl[4][4];
+//   
+//   fak=0.5f;	//TAU moras izracunati le na LMARK1 markerjih
+// 
+//   lp=l+1; lm=l-1;	
+//   if(!(LMARK[lm]&(LMARK1))) {lm=l,fak=1.f;} else if(!(LMARK[lp]&LMARK1)) {lp=l;fak=1.f;};
+//   izrAStressL(lp,taup);									//calculate stress in lxp nad lxm points
+//   izrAStressL(lm,taum);
+//   DTAUl[1][1]=fak*(taup[1][1]-taum[1][1]);							//calculate derivatives of stress point: eg DTAU[2][1]=DTau_yx /dx
+//   DTAUl[2][1]=fak*(taup[2][1]-taum[2][1]);
+//   DTAUl[3][1]=fak*(taup[3][1]-taum[3][1]);
+// 
+//   lp=l+I; lm=l-I;
+//   if(!(LMARK[lm]&LMARK1)) {lm=l,fak=1.f;} else if(!(LMARK[lp]&LMARK1)) {lp=l;fak=1.f;};
+//   izrAStressL(lp,taup);									//calculate stress in lyp nad lym points
+//   izrAStressL(lm,taum);
+//   DTAUl[1][2]=fak*(taup[1][2]-taum[1][2]);							//calculate derivatives of stress point: eg DTAU[2][1]=DTau_yx /dx
+//   DTAUl[2][2]=fak*(taup[2][2]-taum[2][2]);
+//   DTAUl[3][2]=fak*(taup[3][2]-taum[3][2]);
+// 
+//   lp=l+I*J; lm=l-I*J;
+//   if(!(LMARK[lm]&LMARK1)) {lm=l,fak=1.f;} else if(!(LMARK[lp]&LMARK1)) {lp=l;fak=1.f;};
+//   izrAStressL(lp,taup);									//calculate stress in lyp nad lym points
+//   izrAStressL(lm,taum);
+//   DTAUl[1][3]=fak*(taup[1][3]-taum[1][3]);							//calculate derivatives of stress point: eg DTAU[2][1]=DTau_yx /dx
+//   DTAUl[2][3]=fak*(taup[2][3]-taum[2][3]);
+//   DTAUl[3][3]=fak*(taup[3][3]-taum[3][3]);
+// 
+//   divTaul[1]=DTAUl[1][1]+DTAUl[1][2]+DTAUl[1][3];
+//   divTaul[2]=DTAUl[2][1]+DTAUl[2][2]+DTAUl[2][3];
+//   divTaul[3]=DTAUl[3][1]+DTAUl[3][2]+DTAUl[3][3];
+// }
+
+
+//Computes equilibrium distribution function
+void izrFeq(float **feq){
+  float ue,u2;
+  int m;
+  
+  float Ul2[5]; //Local matrices
+
+  #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500) private(Ul2,ue,u2,m)			//Parallelisation -- (static, 20000) is recommended by David Seč (not necessarily the best option)
+  for(int l=0;l<NMAX;l++){
+    if(!(LMARK[l]==LMARK1)) continue;
+    fillUl2(l,Ul2);			//Save global values from U to a local variable Ul2
+    u2=fsqr(Ul2[1]) + fsqr(Ul2[2]) + fsqr(Ul2[3]);	//Velocity squared
+    
+    for(m=0;m<ST_HITROSTI;m++){
+      ue=Ul2[1]*E[m][1] + Ul2[2]*E[m][2] + Ul2[3]*E[m][3];		//Product of velocity and characteristic vectors
+      feq[m][l]=WKONST[m]*Ul2[4]*(1.f+3.f*ue+9.f/2.f*fsqr(ue)-1.5*u2);
+    }
+  }
+}
+
+
+
+//Computes the forcing P
+void izrP(){
+  int l,m;
+  
+  float Ul2[5];
+  float Dxx2[7],Dyy2[7],Dzz2[7],Dx2[7],Dy2[7],Dz2[7]; //Vectors, where the derivatives are stored
+  float divTaul[4];  //Divergence of the stress tensor
+  
+  compute_stress_tensor(); 
+  
+  #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500) private(m,Dxx2,Dyy2,Dzz2,Dx2,Dy2,Dz2,Ul2, divTaul)
+  for(l=0;l<NMAX;l++){
+    if(!(LMARK[l]==LMARK1)) continue;
+    fillUl2(l,Ul2);		//Fill local values of U
+    izrOdvQ2(l,Dxx2,Dyy2,Dzz2,Dx2,Dy2,Dz2);	//Calculate derivatives in Q
+        
+    izrOdvAStressL(l,divTaul);		//Calculate the divergence of the stress tensor
+
+    
+    
+    float Fx=0,Fy=0,Fz=0;
+    for(m=1;m<4;m++){ Fx-=H[m][l]*Dx2[m]; Fy-=H[m][l]*Dy2[m]; Fz-=H[m][l]*Dz2[m];}	//Alternative description of stress tensor (Lisbon version) -- truly incompressible nematic
+    for(m=4;m<7;m++){ Fx-=2*H[m][l]*Dx2[m]; Fy-=2*H[m][l]*Dy2[m]; Fz-=2*H[m][l]*Dz2[m];}
+//     for(m=1;m<4;m++){ Fx-=(Dxx2[m]+Dyy2[m]+Dzz2[m])*Dx2[m]; Fy-=(Dxx2[m]+Dyy2[m]+Dzz2[m])*Dy2[m]; Fz-=(Dxx2[m]+Dyy2[m]+Dzz2[m])*Dz2[m];}
+//     for(m=4;m<7;m++){ Fx-=2*(Dxx2[m]+Dyy2[m]+Dzz2[m])*Dx2[m]; Fy-=2*(Dxx2[m]+Dyy2[m]+Dzz2[m])*Dy2[m]; Fz-=2*(Dxx2[m]+Dyy2[m]+Dzz2[m])*Dz2[m];}
+    float forceX, forceY, forceZ;
+    forceX=divTaul[1] +Fx; FORCE[1][l]=forceX;	//First moment of forcing has contributions from TAU stress tensor and free energy derivatives as well. Velocity is corrected by the factor of FORCE
+    forceY=divTaul[2] +Fy; FORCE[2][l]=forceY;
+    forceZ=divTaul[3] +Fz; FORCE[3][l]=forceZ;
+    for(m=0;m<ST_HITROSTI;m++){
+      float ue=Ul2[1]*E[m][1] + Ul2[2]*E[m][2] + Ul2[3]*E[m][3];		//Product of the velocity and the characteristic vectors
+      float eF=E[m][1]*forceX + E[m][2]*forceY +E[m][3]*forceZ;
+      float uF=Ul2[1]*forceX  + Ul2[2]*forceY  + Ul2[3]*forceZ;
+      P[m][l]=(1.f-DT/2.f/TAUF)*WKONST[m]*( 3.f*eF - 3.f*uF + 9.f *ue*eF );
+      //P[m][l]=0;
+    }
+
+  }
+}
+
+
+
+
+//Transcribes velocity field to a local variable
+void fillUl2(int l, float* ul){
+  ul[1]=U[1][l];
+  ul[2]=U[2][l];
+  ul[3]=U[3][l];
+  ul[4]=U[4][l];
+
+}
+
+
+//Computes streaming location
+int calcLBlnew(int l, int m){
+  int i,j,k;
+  
+  i=i_vr(l)+E[m][1];
+  j=j_vr(l)+E[m][2];
+  k=k_vr(l)+E[m][3];
+
+  return i+j*I+k*I*J;
+}
+
+
+//Computes bounce-back location
+ int calcLBmnew(int m){
+  if(m%2==0) return m-1;		//Opposite vectors are arranged in odd-even order
+  else return m+1;
+}
+
+/*
+//No-slip bounce-back for horizontal plane
+ int calcLBmnew_noslipZ(int m){
+  if(m==5) return 6; else if(m==6) return 5; else if(m==11) return 14; else if(m==14) return 11; else if(m==12) return 13; else if(m==13) return 12; else if(m==15) return 18; else if(m==18) return 15; else if(m==16) return 17; else if(m==17) return 16; else return 0;
+}*/
+
+// Periodic BC
+int calcLpbc(int l){
+  int xp2=0, yp2=0, zp2=0;	//Normal pbc
+   
+  if(i_vr(l)==0) xp2=I-2;
+  if(i_vr(l)==I-1) xp2=-(I-2);
+  if(j_vr(l)==0) yp2=I*(J-2);
+  if(j_vr(l)==J-1) yp2=-I*(J-2);
+  if(k_vr(l)==0) zp2=I*J*(K-2);
+  if(k_vr(l)==K-1) zp2=-I*J*(K-2);
+      
+  return (l+ xp2 + yp2 + zp2);
+}
+
+
+// Open BC
+// int calcLpbc(int l){
+//    int xp2=0, yp2=0, zp2=0;    //Open pbc
+// 
+//    if(i_vr(l)==0) xp2=1;
+//    if(i_vr(l)==I-1) xp2=-1;
+//    if(j_vr(l)==0) yp2=I;
+//    if(j_vr(l)==J-1) yp2=-I;
+//    if(k_vr(l)==0) zp2=I*J;
+//    if(k_vr(l)==K-1) zp2=-I*J;
+// 
+//    return (l+ xp2 + yp2 + zp2);
+// }
+
+//Makes periodic boundary conditions and forces the pressure gradient
+void pbcLB(){
+  #ifdef DELTAP
+  float rho_outletj=0.f, rho_inletj=0.f;
+  int j=0;
+  float rho_outleti=0.f, rho_inleti=0.f;
+  int i=0;
+  for(int l=0;l<NMAX;l++){
+    if(j_vr(l)==1 && (LMARK[l]==LMARK1)){ rho_inletj+=U[4][l]; j++;}
+    else if(j_vr(l)==J-2  && (LMARK[l]==LMARK1)) rho_outletj+=U[4][l];
+    if(i_vr(l)==1 && (LMARK[l]==LMARK1)){ rho_inleti+=U[4][l]; i++;}
+    else if(i_vr(l)==I-2  && (LMARK[l]==LMARK1)) rho_outleti+=U[4][l];
+  }
+  rho_inletj/=j;
+  rho_outletj/=j;
+  rho_inleti/=i;
+  rho_outleti/=i;
+  #endif
+  
+  #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500)
+  for(int l=0; l<NMAX; l++){
+    if(LMARK[l]&LMARKRP){
+      int lpbc=calcLpbc(l); //if(LMARK[l]&LMARK1 && !(LMARK[lpbc]&LMARK1)){ printf("Napacna alokacija in RP!!!"); exit(1);}	
+
+      if(LMARK[l]&LMARK1){
+	float feqnew=0;
+	for(int m=0; m<ST_HITROSTI; m++){
+	  #ifdef DELTAP
+	  if(j_vr(l)==0) feqnew=WKONST[m]*(DENSITYINIT+DELTAPY*(1.f/2.f+1.f/J/2)-rho_outletj);
+	  else if(j_vr(l)==J-1) feqnew=WKONST[m]*(DENSITYINIT-DELTAPY*(1.f/2.f+1.f/J/2)-rho_inletj);
+	  else if(i_vr(l)==0) feqnew=WKONST[m]*(DENSITYINIT+DELTAPX*(1.f/2.f+1.f/I/2)-rho_outleti);
+	  else if(i_vr(l)==I-1) feqnew=WKONST[m]*(DENSITYINIT-DELTAPX*(1.f/2.f+1.f/I/2)-rho_inleti);
+	  #endif
+
+	  P[m][l]=P[m][lpbc]; FEQ[m][l]=FEQ[m][lpbc]+feqnew; F[m][l]=F[m][lpbc]+feqnew;}
+      }
+    }
+  }
+}
+
+//Time step in LB
+void izrCasKorU(){
+  int l,m;
+  int lnew,mnew;		//Where the function is streamed
+  
+  izrFeq(FEQ);		//Computes equilibrium distribution
+  izrP();           //Computes forcing
+  
+  pbcLB();		//PBC
+
+  #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500) private(lnew,m)
+  for(l=0;l<NMAX;l++){
+    if( ((LMARK[l]&LMARK1))){
+      for(m=0;m<ST_HITROSTI;m++){
+	lnew=calcLBlnew(l,m);
+	if(lnew>=0 && lnew<NMAX) FNEW[m][lnew]=F[m][l]-DT/TAUF*(F[m][l]-FEQ[m][l]) + P[m][l];		//Streaming and collision
+      }
+    }
+  }
+  
+  
+  //bounce back -- implementation of surfaces
+  #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500) private(lnew,m,mnew)
+  for(l=0;l<NMAX;l++){	//bounce-back
+    if(!(LMARK[l]&LMARK1) && !(LMARK[l]&LMARK0) ){
+      for(m=1;m<ST_HITROSTI;m++){
+	//mnew=calcLBmnew_noslipZ(m);		//Bounce back at the boundary
+	//lnew=calcLpbc(calcLBlnew(l,mnew));	//To which point the function is returned
+	mnew=calcLBmnew(m);		//Bounce back at the boundary
+	lnew=calcLBlnew(l,mnew);	//To which point the function is returned
+	int lnewE=calcLBlnew(lnew,mnew); //Point E for interpolation
+	if(lnew>=0 && lnew<NMAX && (LMARK[lnew]==(LMARK1)) ) {	//If function is streamed to the bulk
+	  //Bounce back is based on Bouzidi, Phys Fluid 2001	Factor is conneced to the form of FEQ (B)
+//  	  FNEW[mnew][lnew]=FNEW[m][l]-2.f*3.f*WKONST[m]*U[4][lnew]*(E[m][1]*UNIC[1][l]+E[m][2]*UNIC[2][l]+E[m][3]*UNIC[3][l]);		//No interpolation, surface at 1/2
+	  FNEW[mnew][lnew]=0.5f*FNEW[m][l]+0.5f*FNEW[mnew][calcLBlnew(lnew,mnew)]/*;//*/ - 3.f*WKONST[m]*U[4][lnew]*(E[m][1]*UNIC[1][l]+E[m][2]*UNIC[2][l]+E[m][3]*UNIC[3][l]);	//Interpolation at q=1
+	  //FNEW[mnew][lnew]=FNEW[m][lnew];			//Interpolation at q=0
+	}
+      }
+    }
+  }
+ 
+  
+  calcFNEW2F();
+  
+  //Calculate velocity profile
+  #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500)
+  for(l=0;l<NMAX;l++) {calcF2Ul(l,F); /*printf("%d; %f ",LMARK[l], U[1][l]);*/}; 
+  
+  
+  #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500)
+  for(l=0;l<NMAX;l++) for(int m=1; m<5; m++) U[m][l]=U[m][calcLpbc(l)]; 
+}
+
+//The main function for the coupled dynamics
+void izrQU(){
+  int t;
+  char buf[200]; 
+  double fe;
+  FILE *free_energy;
+  float x,y,z,r;
+  
+  //struct ALL ACTIVE DROPLET VARIABLES
+  ActiveDropletVariables act_var;
+  ActiveDropletVariables* active_var = &act_var;
+  active_droplet_var(active_var);
+
+ 
+  //struct init SPHERE
+  ActiveDroplet sphere1;
+  ActiveDroplet* sphere1_ptr = &sphere1;
+  
+  double center1[3] = {0, 0, 0};
+  sphere1.center = center1;
+  sphere1.radius = R;
+  init_droplet(sphere1_ptr);
+  printf("Len after init: %d\n", sphere1.num_points);
+  //     gradU=matrix(1,3,0,sphere1.num_points-1);
+
+
+  //struct FIBONACCI
+  FibonacciLattice sphere_lattice;
+  FibonacciLattice* init_sphere_lattice = &sphere_lattice;
+  fibonacci_init(init_sphere_lattice, 100, center1, R); 
+  printf("Fibonacci num_points: %d\n", sphere_lattice.num_points);
+    
+    
+
+    
+    
+  //Initial time-steps in the Q-tensor. Mostly not needed.
+  for(t=0;t<STCASKORINIT;t++){
+    if(t%STKONV==0){
+    #if TWRITEDIR
+      //sprintf(buf,"./data/%s_dir_%07d-INIT.raw", imeIzpis,t);	zapisQ2DIRraw(buf);
+      sprintf(buf,"./data/%s_op-INIT_%07d.raw",imeIzpis, t);		zapisQ2OPraw(buf);
+    #endif 
+    }
+    izrCasKorQ();	//A time step of the Q-tensor dynamics
+  }
+  
+  
+  //The main dynamics
+  for(t=0;t<STCASKOR;t++){
+     
+    if(t%STKONV==0){      
+    #if TWRITEU
+      sprintf(buf,"./data/%s_u_%07d.vtk", imeIzpis,t);		zapisUraw(buf);
+    #endif
+//     #if TWRITEDEN
+//       sprintf(buf,"./data/%s_den_%07d.vtk", imeIzpis,t);	zapisDENraw(buf);
+//     #endif
+//     #if TWRITEDIR
+//       sprintf(buf,"./data/%s_dir_%07d.vtk", imeIzpis,t);	zapisQ2DIRraw(buf);
+//     #endif 
+    #if TWRITEOP
+      sprintf(buf,"./data/%s_op_%07d.vtk",imeIzpis, t);		zapisQ2OPraw(buf);
+    #endif
+    #if TWRITEQ
+      sprintf(buf,"./data/%s_Q_%07d.raw",imeIzpis, t);		zapisQ(buf);
+    #endif
+    #if TWRITELOC
+      sprintf(buf,"./data/%s_locDefect_%07d.txt",imeIzpis, t);		locDefect(&sphere1,buf);
+      
+      sprintf(buf,"./data/%s__%07d.txt",imeIzpis, t);		locDefect(&sphere1,buf);
+    #endif
+    
+    }
+    
+    if(t%STKONVV==0){      
+    #if TWRITEDEN
+      sprintf(buf,"./data/%s_den_%07d.vtk", imeIzpis,t);	zapisDENraw(buf);
+    #endif
+    #if TWRITEDIR
+      sprintf(buf,"./data/%s_dir_%07d.vtk", imeIzpis,t);	zapisQ2DIRraw(buf);
+    #endif 
+    }
+    
+    if(t%STKONG==0){
+    #if TWRITEgrad
+        zapis_all_grad(&sphere1);
+        zapis_fibonacci_average_velocity(&sphere1, &sphere_lattice, &act_var);
+        zapisDistance();
+    #endif
+    }
+        
+    if(t%STDVF==0){
+    #if TWRITEDVF
+        zapisDVF(buf);
+    #endif   
+    }
+    
+    if(t%STFE==0){
+    #if TWRITEFE
+//         zapisFE(buf);
+        fe=calculate_free_energy();
+        free_energy=fopen("./data/free_energy.txt","a");
+        
+        if(free_energy==NULL){
+            fprintf(stderr, "Failed opening file. a mode. file=%s\n", "./free_energy.txt");
+            exit(1); 
+        }
+        
+        fprintf(free_energy, "%d;%lf\n",t,fe);
+        printf("%d,%lf\n", t,fe);
+        fclose(free_energy);
+    #endif
+    }
+    
+    if(t%STCOPY==0){
+    #if TWRITECOPY
+        sprintf(buf, "%s_DATA.bin", imeIzpis); zapisDATA(buf);
+    #endif    
+    } 
+    
+
+    
+    
+    for(int i=0; i<STCASKORQ; i++) izrCasKorQ(); //Q-tensor time step
+
+
+    
+//hitrostno polje
+    izrCasKorU();	//LB time step
+    
+    
+
+    //surface velocity vX, vY, vZ
+    velocity_sphere(&sphere1, &act_var);
+//     velocity_sphere(&sphere1,0.0001,2);
+
+    velocity_grad(&sphere1);
+//     defect_grad(&sphere1);
+    
+    if(t%STKONV==0){ 
+    #if TWRITEU
+      sprintf(buf,"./data/%s_grad_%07d.vtk", imeIzpis,t);		zapisUgrad(buf);
+//       sprintf(buf,"./data/%s_gradDefect_%07d.vtk", imeIzpis,t);		zapisUgradDEFECT(buf);
+    #endif
+    }
+    
+    
+
+
+
+    
+    
+  }
+}

dipol/fibonacci/saturn_ring/state_A1_test/acn_funk_op.c

@@ -0,0 +1,621 @@
+//Sets up the initial condition for the Q-tensor, QNIC tensor and logical markers LMARK in each datapoint
+void zacPriblQ(long seed){
+  int l;
+  float x,y,z,r;
+  int t;
+ 
+  for(l=0;l<NMAX;l++){
+    x=i_vr(l)-(I-1.f)/2.f;
+    y=j_vr(l)-(J-1.f)/2.f;
+    z=k_vr(l)-(K-1.f)/2.f;
+        
+    r=sqrt(fsqr(x) + fsqr(y) + fsqr(z));
+
+    if(fsqr(x) + fsqr(y) + fsqr(z)<=fsqr(R)) dir2QRP(l,x/r,y/r,z/r,  x/r,y/r,z/r,  0.f,0.f,0.f,LMARKPOV, seed);
+    else if(k_vr(l) == 0 || k_vr(l) == 1) dir2QRP(l,1.f,0.f,0.f,  1.f,0.f,0.f,  0.f,0.f,0.f,LMARKPOV, seed);
+    else if(k_vr(l) == K-1 || k_vr(l) == K-2) dir2QRP(l,-1.f,0.f,0.f,  -1.f,0.f,0.f,  0.f,0.f,0.f,LMARKPOV, seed);
+       
+    else 
+    dir2QBulk(l, seed);
+  }
+  
+  
+  for(l=0;l<NMAX;l++){
+	int lxp,lxm,lyp,lym,lzp,lzm;
+	if(i_vr(l)==0) {lxp=l+1; lxm=l+I-1;}		//periodicnost po x smeri
+	else if(i_vr(l)==I-1) {lxp=l-(I-1); lxm=l-1;}
+	else {lxp=l+1; lxm=l-1;};
+
+	if(j_vr(l)==0) {lyp=l+I; lym=l+(J-1)*I;}	//periodicnost po y smeri
+	else if(j_vr(l)==J-1) {lyp=l-(J-1)*I; lym=l-I;}
+	else {lyp=l+I; lym=l-I;};
+
+	if(k_vr(l)==0) {lzp=l+I*J; lzm=l+(K-1)*I*J;}	//periodicnost po z smeri
+	else if(k_vr(l)==K-1) {lzp=l-(K-1)*I*J; lzm=l-I*J;}
+	else {lzp=l+I*J; lzm=l-I*J;}; 
+    
+	
+	if(LMARK[lxp]&(LMARKPOV|LMARKKOL) && LMARK[lxm]&(LMARKPOV|LMARKKOL) && LMARK[lyp]&(LMARKPOV|LMARKKOL) && LMARK[lym]&(LMARKPOV|LMARKKOL) && LMARK[lzp]&(LMARKPOV|LMARKKOL) && LMARK[lzm]&(LMARKPOV|LMARKKOL)){ LMARK[l]=(LMARKKOL|LMARK0);}
+  }
+
+  for(l=0;l<NMAX;l++){
+    if(i_vr(l)==0 || i_vr(l)==I-1 || j_vr(l)==0 || j_vr(l)==J-1 || k_vr(l)==0 || k_vr(l)==K-1) LMARK[l]=(LMARK[l]|LMARKRP);														
+    else if(LMARK[l]&LMARKKOL) for(int m=1; m<ST_HITROSTI; m++) if(LMARK[calcLBlnew(l, m)]&(LMARK1|LMARKPOV)) LMARK[l]=LMARKKOL;
+  }
+
+
+}
+
+
+//Initialise boundary condition on the surfaces
+void dir2QRP(int l, float dirX, float dirY, float dirZ, float nX, float nY, float nZ, float vX, float vY, float vZ, char lmark, long seed){
+  float seq=SEQ;
+  float x,y,z;
+  int count=0;
+  
+  LMARK[l]=lmark;
+  
+  UNIC[1][l]=vX;		//surface velocity
+  UNIC[2][l]=vY;
+  UNIC[3][l]=vZ;
+  
+ 
+  float norma=sqrt(fsqr(nX)+fsqr(nY)+fsqr(nZ));
+  if(norma==0) norma=1.f;
+  QNIC[4][l]=nX/norma;
+  QNIC[5][l]=nY/norma;
+  QNIC[6][l]=nZ/norma;
+  
+
+  x=(float)(i_vr(l))-(float)(I-1)/2+0.001;
+  y=(float)(j_vr(l))-(float)(J-1)/2+0.001;
+  z=(float)(k_vr(l))-(float)(K-1)/2+0.001;
+//   if(fsqr(x)+fsqr(y)<fsqr(I/2-3)) {dirZ=cos(z/K*M_PI),dirY=0,dirX=sin(z/K*M_PI);}
+//   else {dirX=1; dirY=0; dirZ=0;};
+  
+//   norma=sqrt(fsqr(dirX)+fsqr(dirY)+fsqr(dirZ));
+//   if(norma==0) norma=1.f;
+//   dirX=dirX/norma;//x/sqrt(x*x+y*y+z*z);
+//   dirY=dirY/norma;//y/sqrt(x*x+y*y+z*z);
+//   dirZ=dirZ/norma;//z/sqrt(x*x+y*y+z*z);
+//   
+
+//   vX=y/sqrt(x*x+y*y+z*z);
+  
+  
+  
+//   dirX=x/sqrt(x*x+y*y+z*z);
+//   dirY=y/sqrt(x*x+y*y+z*z);
+//   dirZ=z/sqrt(x*x+y*y+z*z);
+  
+
+  QNIC[1][l]=dirX;
+  QNIC[2][l]=dirY;
+  QNIC[3][l]=dirZ;
+//   dirX=0; dirY=0; dirZ=1;
+  
+  
+  
+//   //random director
+//   float phi, theta;
+//   phi=ran2(&seed)*2*M_PI;
+//   theta=acos(2*(ran2(&seed))-1);
+//   
+// //   x=cos(phi)*sin(theta);
+// //   y=sin(phi)*sin(theta);
+// //   z=cos(theta);
+//   
+//   dirX=x/sqrt(x*x+y*y+z*z);
+//   dirY=y/sqrt(x*x+y*y+z*z);
+//   dirZ=z/sqrt(x*x+y*y+z*z);
+  
+ 
+  
+  Q[1][l]=0.5f*seq*(3*dirX*dirX-1);
+  Q[2][l]=0.5f*seq*(3*dirY*dirY-1);
+  Q[3][l]=0.5f*seq*(3*dirZ*dirZ-1);
+  Q[4][l]=1.5f*seq*(dirX*dirY);
+  Q[5][l]=1.5f*seq*(dirX*dirZ);
+  Q[6][l]=1.5f*seq*(dirY*dirZ);
+}
+
+
+//Initialise boundary condition in the bulk
+void dir2QBulk(int l,long seed){
+  LMARK[l]=LMARK1;
+  float dirX, dirY, dirZ;
+  float seq=SEQ;//ran2(&seed)*1e-30;
+  float x,y,z;
+  
+  
+//   x=(float)(i_vr(l))-(float)(I-1)/2+0.001;
+//   y=(float)(j_vr(l))-(float)(J-1)/2+0.001;
+//   z=(float)(k_vr(l))-(float)(K-1)/2+0.001;
+
+//   //loop
+//   if(fsqr(x)+fsqr(y)<fsqr(1/8)) {dirZ=cos(z/K*M_PI),dirY=0,dirX=sin(z/K*M_PI);}
+//   dirX=0; dirY=0; dirZ=1;
+  
+
+  float q, e_x, r;
+  
+  q=1.2;
+  e_x=1.0;
+  
+  
+  
+  
+  x=(float)(i_vr(l))-(float)(I-1)/2+0.001;
+  y=(float)(j_vr(l))-(float)(J-1)/2+0.001;
+  z=(float)(k_vr(l))-(float)(K-1)/2+0.001;
+  
+  r=sqrt(x*x+y*y+z*z)*sqrt(x*x+y*y+z*z)*sqrt(x*x+y*y+z*z);
+  
+  dirX=e_x+q*R*R*x/r;
+  dirY=+q*R*R*y/r;
+  dirZ=+q*R*R*z/r;
+ 
+  
+//   x=1;
+//   y=0;
+//   z=1;
+//   
+//   dirX=x/sqrt(x*x+y*y+z*z);
+//   dirY=y/sqrt(x*x+y*y+z*z);
+//   dirZ=z/sqrt(x*x+y*y+z*z);
+  
+  
+//   
+// 
+//   QNIC[1][l]=dirX;
+//   QNIC[2][l]=dirY;
+//   QNIC[3][l]=dirZ;
+  
+  
+  Q[1][l]=0.5f*seq*(3*dirX*dirX-1);
+  Q[2][l]=0.5f*seq*(3*dirY*dirY-1);
+  Q[3][l]=0.5f*seq*(3*dirZ*dirZ-1);
+  Q[4][l]=1.5f*seq*(dirX*dirY);
+  Q[5][l]=1.5f*seq*(dirX*dirZ);
+  Q[6][l]=1.5f*seq*(dirY*dirZ);
+  
+  UNIC[1][l]=0;		//in the bulk, surface velocity is zero
+  UNIC[2][l]=0;
+  UNIC[3][l]=0;
+
+}
+
+
+//Compute local Q-tensor derivatives
+void izrOdvQ2(int l, float* dxx, float* dyy, float* dzz, float* dx, float* dy, float* dz){
+  int i,lxp,lxm,lyp,lym,lzp,lzm;
+
+	
+	lxp=l+1; lxm=l-1;
+	lyp=l+I; lym=l-I;
+	lzp=l+I*J; lzm=l-I*J;
+
+
+	for(i=1;i<=6;i++)				//derivatives of Q. Watch out: 2. derivatives are calculated only in the bulk
+	{
+	  
+		dxx[i]=(Q[i][lxp]+Q[i][lxm]-2*Q[i][l]);
+		dyy[i]=(Q[i][lyp]+Q[i][lym]-2*Q[i][l]);
+		dzz[i]=(Q[i][lzp]+Q[i][lzm]-2*Q[i][l]);
+	  
+
+		dx[i]=0.5f*(Q[i][lxp]-Q[i][lxm]);
+		dy[i]=0.5f*(Q[i][lyp]-Q[i][lym]);
+		dz[i]=0.5f*(Q[i][lzp]-Q[i][lzm]);
+	};
+
+	
+}
+
+//Compute forward derivative of the Q-tensor on the surface
+void izrOdvQ2POV(int l, float* dx, float* dy, float* dz){
+  int i,lxp,lxm,lyp,lym,lzp,lzm,lsxp,lsxm,lsyp,lsym,lszp,lszm;
+	float fakX=0.5,fakY=0.5,fakZ=0.5;
+	
+	lxp=l+1; lxm=l-1;
+	lyp=l+I; lym=l-I;
+	lzp=l+I*J; lzm=l-I*J;
+
+	lsxp=lxp; lsxm=lxm; lsyp=lyp; lsym=lym; lszp=lzp; lszm=lzm;
+
+
+		if(LMARK[lxm]&LMARKKOL) {lsxp=lxp;lsxm=l;fakX=1.;}
+		else if(LMARK[lxp]&LMARKKOL) {lsxp=l;lsxm=lxm;fakX=1.;};
+
+		if(LMARK[lym]&LMARKKOL) {lsyp=lyp;lsym=l;fakY=1.;}
+		else if(LMARK[lyp]&LMARKKOL) {lsyp=l;lsym=lym;fakY=1.;};
+
+		if(LMARK[lzm]&LMARKKOL) {lszp=lzp;lszm=l;fakZ=1.;}
+		else if(LMARK[lzp]&LMARKKOL) {lszp=l;lszm=lzm;fakZ=1.;};
+
+
+	for(i=1;i<=6;i++)				//derivatives of Q
+	{
+
+		dx[i]=fakX*(Q[i][lsxp]-Q[i][lsxm]);
+		dy[i]=fakY*(Q[i][lsyp]-Q[i][lsym]);
+		dz[i]=fakZ*(Q[i][lszp]-Q[i][lszm]);
+	};
+
+	
+}
+
+//free energy
+double calculate_free_energy(){
+  
+  double free_energy=0;
+  #pragma omp parallel for reduction (+:free_energy)
+  for(int l=0;l<NMAX;l++){
+      
+    if(LMARK[l]==LMARK1){
+      float Dxx2[7],Dyy2[7],Dzz2[7],Dx2[7],Dy2[7],Dz2[7],Dxy2[7],Dxz2[7],Dyz2[7];
+      izrOdvQ2(l,Dxx2,Dyy2,Dzz2,Dx2,Dy2,Dz2);
+//       izrOdvQ2(l,Dxx2,Dyy2,Dzz2,Dxy2,Dxz2,Dyz2,Dx2,Dy2,Dz2);
+
+      free_energy+=
+	//L1 part
+	0.5*( sqr(Dx2[1]) + sqr(Dx2[2]) +sqr(Dx2[3]) + 2*sqr(Dx2[4]) + 2*sqr(Dx2[5]) + 2*sqr(Dx2[6]) 
+	+ sqr(Dy2[1]) + sqr(Dy2[2]) +sqr(Dy2[3]) + 2*sqr(Dy2[4]) + 2*sqr(Dy2[5]) + 2*sqr(Dy2[6])
+        + sqr(Dz2[1]) + sqr(Dz2[2]) +sqr(Dz2[3]) + 2*sqr(Dz2[4]) + 2*sqr(Dz2[5]) + 2*sqr(Dz2[6]))
+// 	//L2 part
+// 	+0.5*L21*(sqr(Dx2[1]) +sqr(Dx2[4]) +sqr(Dx2[5]) +sqr(Dy2[2]) +sqr(Dy2[4]) +sqr(Dy2[6]) +sqr(Dz2[3]) +sqr(Dz2[5]) +sqr(Dz2[6]) +2*Dz2[6]*Dy2[2] +2*Dz2[5]*Dy2[3] +2*Dz2[3]*Dy2[6] +2*Dz2[5]*Dx2[1] +2*Dy2[4]*Dx2[1] +2*Dz2[6]*Dx2[4] +2*Dy2[2]*Dx2[4] +2*Dz2[3]*Dx2[5] +2*Dy2[6]*Dx2[5]	      )
+// 	//Cholesteric part
+// 	+2*q0*(
+// 		+Q[4][l]*Dz2[1]
+// 		-Q[4][l]*Dz2[2]
+// 		+Q[2][l]*Dz2[4]
+// 		-Q[1][l]*Dz2[4]
+// 		+Q[6][l]*Dz2[5]
+// 		-Q[5][l]*Dz2[6]
+// 		
+// 		+Q[5][l]*Dy2[3]
+// 		-Q[5][l]*Dy2[1]
+// 		+Q[1][l]*Dy2[5]
+// 		-Q[6][l]*Dy2[4]
+// 		+Q[4][l]*Dy2[6]
+// 		-Q[3][l]*Dy2[5]
+// 
+// 		+Q[6][l]*Dx2[2]
+// 		-Q[6][l]*Dx2[3]
+// 		+Q[5][l]*Dx2[4]
+// 		-Q[4][l]*Dx2[5]
+// 		+Q[3][l]*Dx2[6]
+// 		-Q[2][l]*Dx2[6]
+// 
+// 		)
+	//Phase part
+        + 0.5*A*(sqr(Q[1][l])+sqr(Q[2][l])+sqr(Q[3][l])+2*sqr(Q[4][l])+2*sqr(Q[5][l])+2*sqr(Q[6][l]))
+					+B/3.*(Q[1][l]*(sqr(Q[1][l])+sqr(Q[4][l])+sqr(Q[5][l]))+Q[2][l]*(sqr(Q[2][l])+sqr(Q[4][l])+sqr(Q[6][l]))+Q[3][l]*(sqr(Q[3][l])+sqr(Q[5][l])+sqr(Q[6][l]))+2*Q[6][l]*(Q[4][l]*Q[5][l]+Q[2][l]*Q[6][l]+Q[3][l]*Q[6][l])+2*Q[5][l]*(Q[1][l]*Q[5][l]+Q[3][l]*Q[5][l]+Q[4][l]*Q[6][l])+2*Q[4][l]*(Q[1][l]*Q[4][l]+Q[2][l]*Q[4][l]+Q[5][l]*Q[6][l]))
+					+0.25*C*sqr(sqr(Q[1][l])+sqr(Q[2][l])+sqr(Q[3][l])+2*sqr(Q[4][l])+2*sqr(Q[5][l])+2*sqr(Q[6][l]))
+					-(0.75*A*sqr(SEQ)+0.25*B*SEQ*sqr(SEQ)+9*C*sqr(sqr(SEQ))/16.0)	
+      ;
+    }
+  }
+//   printf("\nthe free energy: %.20lf",free_energy);
+  return free_energy;
+  
+}
+//
+
+//Computes all needed derivatives of the velocity field in a given point. The result is in units of RESP/TAU and not in the units of LB code (the difference is DT)
+void izrOdvU2(int l, float *ux, float *uy, float *uz){
+	int i,lxp,lxm,lyp,lym,lzp,lzm;
+	
+	float fakX=.5f, fakY=.5f, fakZ=.5f;
+	
+	lxp=l+1; lxm=l-1;
+	lyp=l+I; lym=l-I;
+	lzp=l+I*J; lzm=l-I*J;
+	
+	if(!(LMARK[lxp]&LMARK1)){ lxp=l; fakX=1.f;}	//Forward difference is used because due to TAUF>>DT the approximation of U=0 at the surface is ot fully justified.
+	else if(!(LMARK[lxm]&LMARK1)){ lxm=l; fakX=1.f;}
+	if(!(LMARK[lyp]&LMARK1)){ lyp=l; fakY=1.f;}
+	else if(!(LMARK[lym]&LMARK1)){ lym=l; fakY=1.f;}
+	if(!(LMARK[lzp]&LMARK1)){ lzp=l; fakZ=1.f;}
+	else if(!(LMARK[lzm]&LMARK1)){ lzm=l; fakZ=1.f;}
+	
+	for(i=1;i<=3;i++)				/*derivatives of velocity U*/
+	{
+		ux[i]=fakX*(U[i][lxp]-U[i][lxm])/DT;	/*e.g. Ux[i] is derivative with respect to x of i-th component of U*/
+		uy[i]=fakY*(U[i][lyp]-U[i][lym])/DT;
+		uz[i]=fakZ*(U[i][lzp]-U[i][lzm])/DT;
+	};
+    
+	
+}
+
+int ijk_to_l(int i, int j, int k){
+    return i+j*I+k*I*J;
+}
+
+void interpoliraj_U(int i, int j, int k, float is, float jf, float kf, float* u){
+    for(int m=1;m<=3;m++){
+        u[m]=U[m][ijk_to_l(i,j,k)]*(1-is)*(1-jf)*(1-kf)+
+        U[m][ijk_to_l(i+1,j,k)]*(is)*(1-jf)*(1-kf)+
+        U[m][ijk_to_l(i,j+1,k)]*(1-is)*(jf)*(1-kf)+
+        U[m][ijk_to_l(i+1,j+1,k)]*(is)*(jf)*(1-kf)+
+        U[m][ijk_to_l(i,j,k+1)]*(1-is)*(1-jf)*(kf)+
+        U[m][ijk_to_l(i+1,j,k+1)]*(is)*(1-jf)*(kf)+
+        U[m][ijk_to_l(i,j+1,k+1)]*(1-is)*(jf)*(kf)+
+        U[m][ijk_to_l(i+1,j+1,k+1)]*(is)*(jf)*(kf);
+    }
+}
+
+//Computes all needed derivatives of the velocity field in a given point. The result is in units of RESP/TAU and not in the units of LB code (the difference is DT)
+void izrOdvU2smeri(float* pos, float* n, float *grad){
+    const float korak1=1.1;
+    const float korak2=2.1;
+    
+    float pos1[4];
+    float pos2[4];
+    for(int i=1;i<=3;i++){
+        pos1[i]=pos[i]+n[i]*korak1;
+        pos2[i]=pos[i]+n[i]*korak2;
+    }
+
+    int i1,j1,k1;//indeks točke
+    float if1,jf1,kf1;//ostanki med 0 in 1
+    int i2,j2,k2;//indeks točke + normala*korak
+    float if2,jf2,kf2;
+    
+    i1=(int)pos1[1];
+    j1=(int)pos1[2];
+    k1=(int)pos1[3];
+    
+    if1=pos1[1]-i1;
+    jf1=pos1[2]-j1;
+    kf1=pos1[3]-k1;
+    
+    i2=(int)pos2[1];
+    j2=(int)pos2[2];
+    k2=(int)pos2[3];
+    
+    if2=pos2[1]-i2;
+    jf2=pos2[2]-j2;
+    kf2=pos2[3]-k2;
+    
+    float u1[4];
+    float u2[4];
+    interpoliraj_U(i1,j1,k1,if1,jf1,kf1,u1);
+    interpoliraj_U(i2,j2,k2,if2,jf2,kf2,u2);
+    for(int i=1;i<=3;i++){
+        grad[i]=(u2[i]-u1[i])/(korak2-korak1)/DT;
+    }
+	
+}
+
+
+//Fills in the local Q-tensor from the global matrix
+void fillQl2(int l, float* ql){
+ 	ql[1]=Q[1][l];
+	ql[2]=Q[2][l];
+	ql[3]=Q[3][l];
+	ql[4]=Q[4][l];
+	ql[5]=Q[5][l];
+	ql[6]=Q[6][l];
+}
+
+void fillQlm2(int l, float (*qlm)[4]){
+	qlm[1][1]=Q[1][l];
+	qlm[1][2]=Q[4][l];
+	qlm[1][3]=Q[5][l];
+	qlm[2][1]=Q[4][l];
+	qlm[2][2]=Q[2][l];
+	qlm[2][3]=Q[6][l];
+	qlm[3][1]=Q[5][l];
+	qlm[3][2]=Q[6][l];
+	qlm[3][3]=Q[3][l];
+}
+
+//Computes the molecular field  -- U1 contribution
+void calcU1(int l,float *u1, float* dxx, float* dyy, float* dzz, float* dx, float* dy, float* dz, float* ql){
+  float	tr3,konst=fsqr(ql[1])+fsqr(ql[2])+fsqr(ql[3])+2*(fsqr(ql[4])+fsqr(ql[5])+fsqr(ql[6]));
+
+  u1[4]=dxx[4]+dyy[4]+dzz[4]-2*C*ql[4]*fsqr(ql[4])-B*ql[5]*ql[6]-ql[4]*(A+B*ql[1]+C*fsqr(ql[1])+C*fsqr(ql[2])+B*ql[2]+C*fsqr(ql[3])+2*C*fsqr(ql[5])+2*C*fsqr(ql[6]));
+  u1[5]=dxx[5]+dyy[5]+dzz[5]-2*C*ql[5]*fsqr(ql[5])-B*ql[4]*ql[6]-ql[5]*(A+B*ql[1]+C*fsqr(ql[1])+C*fsqr(ql[2])+B*ql[3]+C*fsqr(ql[3])+2*C*fsqr(ql[4])+2*C*fsqr(ql[6]));
+  u1[6]=dxx[6]+dyy[6]+dzz[6]-2*C*ql[6]*fsqr(ql[6])-B*ql[4]*ql[5]-ql[6]*(A+B*ql[2]+C*fsqr(ql[1])+C*fsqr(ql[2])+B*ql[3]+C*fsqr(ql[3])+2*C*fsqr(ql[4])+2*C*fsqr(ql[5]));
+  
+  u1[1]=dxx[1]+dyy[1]+dzz[1]-A*ql[1]-B*(fsqr(ql[1])+fsqr(ql[4])+fsqr(ql[5]))-C*ql[1]*konst;
+  u1[2]=dxx[2]+dyy[2]+dzz[2]-A*ql[2]-B*(fsqr(ql[2])+fsqr(ql[4])+fsqr(ql[6]))-C*ql[2]*konst;
+  u1[3]=dxx[3]+dyy[3]+dzz[3]-A*ql[3]-B*(fsqr(ql[3])+fsqr(ql[5])+fsqr(ql[6]))-C*ql[3]*konst;
+  
+		tr3=(u1[1]+u1[2]+u1[3])/3.f;	//condition TrQ=0
+		u1[1]-=tr3;
+		u1[2]-=tr3;
+		u1[3]-=tr3;
+
+
+		H[1][l]=u1[1];			//store molecular field which will be needed later for compuing of symmetric and antisymmetric stress tensor
+		H[2][l]=u1[2];
+		H[3][l]=u1[3];
+		H[4][l]=u1[4];
+		H[5][l]=u1[5];
+		H[6][l]=u1[6];
+}
+
+
+
+//Computes U2 contribution to the Q-tensor dynamics -- S_ij
+void calcU2(float *u2, float *ux, float *uy, float *uz, float (*qlm)[4])
+{
+
+  
+	float tr;
+	float D[4][4], O[4][4];
+
+	D[1][1]=ux[1];			//calculate matrix Dij =(Du_i/Dx_j+Du_j/Dx_i)/2  locally 
+	D[1][2]=0.5f*(uy[1]+ux[2]);
+	D[1][3]=0.5f*(uz[1]+ux[3]);
+	D[2][1]=D[1][2];
+	D[2][2]=uy[2];
+	D[2][3]=0.5f*(uz[2]+uy[3]);
+	D[3][1]=D[1][3];
+	D[3][2]=D[2][3];
+	D[3][3]=uz[3];
+
+	O[1][1]=0.f;			//calculate matrix Omega_ij =(Du_i/Dx_j-Du_j/Dx_i)/2  locally
+	O[1][2]=0.5f*(uy[1]-ux[2]);
+	O[1][3]=0.5f*(uz[1]-ux[3]);
+	O[2][1]=-O[1][2];
+	O[2][2]=0.f;
+	O[2][3]=0.5f*(uz[2]-uy[3]);
+	O[3][1]=-O[1][3];
+	O[3][2]=-O[2][3];
+	O[3][3]=0.f;
+
+	tr=	qlm[1][1]*ux[1]+qlm[1][2]*ux[2]+qlm[1][3]*ux[3]+
+		qlm[2][1]*uy[1]+qlm[2][2]*uy[2]+qlm[2][3]*uy[3]+
+		qlm[3][1]*uz[1]+qlm[3][2]*uz[2]+qlm[3][3]*uz[3];
+
+	float D_qlm[4][4];	AMatrixB(D_qlm,D,qlm);
+	float qlm_D[4][4];	AMatrixB(qlm_D,qlm,D);
+	float O_qlm[4][4];	AMatrixB(O_qlm,O,qlm);
+	float qlm_O[4][4]; 	AMatrixB(qlm_O,qlm,O);
+
+
+	u2[1]= XI*(D_qlm[1][1]+qlm_D[1][1]+2.f/3.f*D[1][1])+O_qlm[1][1]-qlm_O[1][1]-2.f*XI*(qlm[1][1]+1.f/3.f)*tr;
+	u2[2]= XI*(D_qlm[2][2]+qlm_D[2][2]+2.f/3.f*D[2][2])+O_qlm[2][2]-qlm_O[2][2]-2.f*XI*(qlm[2][2]+1.f/3.f)*tr;
+	u2[3]= XI*(D_qlm[3][3]+qlm_D[3][3]+2.f/3.f*D[3][3])+O_qlm[3][3]-qlm_O[3][3]-2.f*XI*(qlm[3][3]+1.f/3.f)*tr;
+	u2[4]= XI*(D_qlm[1][2]+qlm_D[1][2]+2.f/3.f*D[1][2])+O_qlm[1][2]-qlm_O[1][2]-2.f*XI*(qlm[1][2])*tr;
+	u2[5]= XI*(D_qlm[1][3]+qlm_D[1][3]+2.f/3.f*D[1][3])+O_qlm[1][3]-qlm_O[1][3]-2.f*XI*(qlm[1][3])*tr;
+	u2[6]= XI*(D_qlm[2][3]+qlm_D[2][3]+2.f/3.f*D[2][3])+O_qlm[2][3]-qlm_O[2][3]-2.f*XI*(qlm[2][3])*tr;
+	
+}
+
+
+//Computes U3 contribution to Q-tensor dynamics -- advection term (v x nabla)Q
+void calcU3(float *u3, float* dx, float* dy, float* dz, float* ul)
+{
+	u3[1]=(-ul[1]*dx[1]-ul[2]*dy[1]-ul[3]*dz[1])/DT;
+	u3[2]=(-ul[1]*dx[2]-ul[2]*dy[2]-ul[3]*dz[2])/DT;
+	u3[3]=(-ul[1]*dx[3]-ul[2]*dy[3]-ul[3]*dz[3])/DT;
+	u3[4]=(-ul[1]*dx[4]-ul[2]*dy[4]-ul[3]*dz[4])/DT;
+	u3[5]=(-ul[1]*dx[5]-ul[2]*dy[5]-ul[3]*dz[5])/DT;
+	u3[6]=(-ul[1]*dx[6]-ul[2]*dy[6]-ul[3]*dz[6])/DT;
+}
+
+
+void izrCasKorQ(){
+
+  
+  #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500)
+  for(int l=0;l<NMAX;l++){
+    if(LMARK[l]==LMARK1){
+      float Dxx2[7],Dyy2[7],Dzz2[7],Dx2[7],Dy2[7],Dz2[7],Ux2[4],Uy2[4],Uz2[4]; //Vectors in which derivatives are stored
+      float Ql2[7];
+      float Ul2[5];
+      float Qlm2[4][4];
+      float u1[7], u2[7], u3[7];
+      fillQl2(l,Ql2);	//Fill the values of Q and U in local matrices
+      fillQlm2(l,Qlm2);
+      fillUl2(l, Ul2);
+      izrOdvQ2(l,Dxx2,Dyy2,Dzz2,Dx2,Dy2,Dz2);	//Computes the derivatives of Q and U
+      izrOdvU2(l,Ux2,Uy2,Uz2);
+      calcU1(l,u1,Dxx2,Dyy2,Dzz2,Dx2,Dy2,Dz2,Ql2);			//calculate u1 contribution -- molecular field H_ij
+      calcU2(u2, Ux2, Uy2, Uz2, Qlm2);					//calculate u1 contribution -- S_ij
+ /*     float Sij=calcU2(u2, Ux2, Uy2, Uz2, Qlm2);
+      printf("Sij = %f",Sij);
+ */     
+      calcU3(u3, Dx2, Dy2, Dz2, Ul2);					//calculate u1 contribution -- advection
+      
+      
+      for(int m=1; m<7; m++) {Ql2[m]+=  DT/STCASKORQ*(u1[m]+u2[m]+u3[m]);}	//Time step in Q
+      float tr3=(Ql2[1]+Ql2[2]+Ql2[3])/3.f;	//Make the trace 0
+      Ql2[1]-=tr3; Ql2[2]-=tr3; Ql2[3]-=tr3; 
+      for(int m=1; m<7; m++) QNEW[m][l]=Ql2[m];	//Write to global matrix
+    }
+    else if(LMARK[l]==LMARKPOV){	//Homeotropic anchoring
+      for(int m=1; m<7; m++) {QNEW[m][l]= Q[m][l];}
+      continue;	//Infinite anchoring
+      float Dx2[7],Dy2[7],Dz2[7]; //Vectors, where the derivatives are stored
+      float n1,n2,n3,qnic[7];
+      float seq=SEQ;
+
+      izrOdvQ2POV(l,Dx2,Dy2,Dz2);
+      
+      
+      n1=QNIC[4][l];	//Surface normal
+      n2=QNIC[5][l];
+      n3=QNIC[6][l];
+      
+      
+            
+//       n1=QNIC[1][l];	//Surface normal
+//       n2=QNIC[2][l];
+//       n3=QNIC[3][l];
+      
+      //Preferred Q-tensor at the surface
+      float dirX, dirY, dirZ;
+      dirX=QNIC[1][l];
+      dirY=QNIC[2][l];
+      dirZ=QNIC[3][l];
+      qnic[1]=0.5f*seq*(3*dirX*dirX-1);
+      qnic[2]=0.5f*seq*(3*dirY*dirY-1);
+      qnic[3]=0.5f*seq*(3*dirZ*dirZ-1);
+      qnic[4]=1.5f*seq*(dirX*dirY);
+      qnic[5]=1.5f*seq*(dirX*dirZ);
+      qnic[6]=1.5f*seq*(dirY*dirZ);
+      
+      
+      for(int m=1; m<7; m++) {QNEW[m][l]= Q[m][l] +  DT/STCASKORQ*(/*2.f**/(n1*Dx2[m]+n2*Dy2[m]+n3*Dz2[m])-WHOM*(Q[m][l]-qnic[m]) ); }	//Time step at the surface
+    
+      float tr3=(QNEW[1][l]+QNEW[2][l]+QNEW[3][l])/3.f;		//Make trace Q zero
+      QNEW[1][l] -= tr3; QNEW[2][l] -= tr3; QNEW[3][l] -= tr3;
+    }
+//     else if(LMARK[l]==LMARKPOV){	//Planar degenerate anchoring
+//       float Dx2[7],Dy2[7],Dz2[7];
+//       float n1,n2,n3,qnic[7];
+// 
+//       izrOdvQ2POV(l,Dx2,Dy2,Dz2);
+//       
+//       
+//       n1=QNIC[4][l];	//Surface normal
+//       n2=QNIC[5][l];
+//       n3=QNIC[6][l];
+//       
+//       
+//       
+//       
+//       float q1=Q[1][l], q2=Q[2][l], q3=Q[3][l], q4=Q[4][l], q5=Q[5][l], q6=Q[6][l]; //local Q
+//       float sq=sqrtf( 2.f/3.f * (fsqr(q1)+fsqr(q2)+fsqr(q3)+2*fsqr(q4)+2*fsqr(q5)+2*fsqr(q6)) ); //local order parameter S
+//       float n11=n1*n1, n22=n2*n2, n33=n3*n3, n12=n1*n2, n13=n1*n3, n23=n2*n3;
+//       float nqn=q1*n11+q2*n22+q3*n33+2*(q4*n12+q5*n13+q6*n23);
+//       float QPOV[7];
+//         QPOV[1]=(sq-2*nqn)*n11   +4*(n11*q1+n12*q4+n13*q5);
+//         QPOV[2]=(sq-2*nqn)*n22   +4*(n12*q4+n22*q2+n23*q6);
+//         QPOV[3]=(sq-2*nqn)*n33   +4*(n13*q5+n23*q6+n33*q3);
+//         QPOV[4]=(sq-2*nqn)*n12 +2*(n11*q4+n12*q2+n13*q6+n12*q1+n22*q4+n23*q5);
+//         QPOV[5]=(sq-2*nqn)*n13 +2*(n11*q5+n12*q6+n13*q3+n13*q1+n23*q4+n33*q5);
+//         QPOV[6]=(sq-2*nqn)*n23 +2*(n12*q5+n22*q6+n23*q3+n13*q4+n23*q2+n33*q6);
+//       
+//       
+//       for(int m=1; m<7; m++) {QNEW[m][l]=Q[m][l]+  DT/STCASKORQ*(2.f*(n1*Dx2[m]+n2*Dy2[m]+n3*Dz2[m])-WHOM*QPOV[m] ); }	//Time step at the surface
+//     
+//       float tr3=(QNEW[1][l]+QNEW[2][l]+QNEW[3][l])/3.f;		//Make trace Q zero
+//       QNEW[1][l] -= tr3; QNEW[2][l] -= tr3; QNEW[3][l] -= tr3;
+//     }
+//     
+    
+    
+  }
+  calcQNEW2Q();		//Set the periodic boundary conditions
+
+  
+  
+  
+}
+
+//Set the periodic boundary conditions
+void calcQNEW2Q(){
+  #pragma omp parallel for num_threads(STPROC) schedule(dynamic, 2500)
+  for(int l=0; l<NMAX; l++){
+   int lpbc=calcLpbc(l);
+   /*if(lpbc!=l) */for(int m=1; m<7; m++) Q[m][l]=QNEW[m][lpbc];
+ }
+}
+

dipol/fibonacci/saturn_ring/state_A1_test/acn_funk_zapis.c

@@ -0,0 +1,637 @@
+
+//Write the necessay data that can be later used for initialisation
+void zapisDATA(char *ime)
+{
+	FILE *file;
+	int l,m;
+	float delovni;
+	char delovni2;
+
+	file=fopen(ime,"w");
+    if(file==NULL){
+        fprintf(stderr, "Failed opening file. w mode. file=%s\n", ime);
+        exit(1);
+    }
+	for(l=0;l<NMAX;l++){
+		for(m=0;m<ST_HITROSTI;m++)
+		{
+			delovni=F[m][l];
+			fwrite(&delovni,sizeof(float),1,file);
+		};
+	}
+	for(l=0;l<NMAX;l++){
+		for(m=1;m<7;m++)
+		{
+			delovni=Q[m][l];
+			fwrite(&delovni,sizeof(float),1,file);
+		};
+	}
+	for(l=0;l<NMAX;l++){
+		for(m=1;m<7;m++)
+		{
+			delovni=QNIC[m][l];
+			fwrite(&delovni,sizeof(float),1,file);
+		};
+	}
+	for(l=0;l<NMAX;l++){
+		{
+			delovni2=LMARK[l];
+			fwrite(&delovni2,sizeof(char),1,file);
+		};
+	}
+	fclose(file);
+}
+
+
+//Write the Q-tensor
+void zapisQ(char *ime)
+{
+	FILE *file;
+	int l,m;
+	float delovni;
+
+	file=fopen(ime,"w");
+    if(file==NULL){
+        fprintf(stderr, "Failed opening file. w mode. file=%s\n", ime);
+         exit(1);
+    }
+
+	for(l=0;l<NMAX;l++){
+		for(m=1;m<7;m++)
+		{
+			delovni=Q[m][l];
+			fwrite(&delovni,sizeof(float),1,file);
+		};
+	}
+
+	fclose(file);
+}
+
+
+//Reads the data file
+void beriDATA(char *ime)
+{
+	FILE *file;
+	int l,m;
+	float delovni;
+	char delovni2;
+	int result;
+
+	file=fopen(ime,"r");
+    if(file==NULL){
+        fprintf(stderr, "Failed opening file. r mode. file=%s\n", ime);
+        exit(1);
+    }
+    
+	for(l=0;l<NMAX;l++){
+		for(m=0;m<ST_HITROSTI;m++)
+		{
+		  result=fread(&delovni,sizeof(float),1,file);
+		  if (result != 1) {fputs ("Error reading data!\n",stderr); exit (2);}
+		  F[m][l]=delovni;
+		};
+	}
+	for(l=0;l<NMAX;l++){
+		for(m=1;m<7;m++)
+		{
+		  result=fread(&delovni,sizeof(float),1,file);
+		  if (result != 1) {fputs ("Error reading data!\n",stderr); exit (2);}
+		  Q[m][l]=delovni;
+		};
+	}
+	for(l=0;l<NMAX;l++){
+		for(m=1;m<7;m++)
+		{
+		  result=fread(&delovni,sizeof(float),1,file);
+		  if (result != 1) {fputs ("Error reading data!\n",stderr); exit (2);}
+		  QNIC[m][l]=delovni;
+		};
+	}
+	for(l=0;l<NMAX;l++){
+		{
+		  result=fread(&delovni2,sizeof(char),1,file);
+		  if (result != 1) {fputs ("Error reading data!\n",stderr); exit (2);}
+		  LMARK[l]=delovni2;
+		}
+	}
+	fclose(file);
+}
+
+
+//Writes the velocity profile along the chosen direction
+void za_profilU(char *ime, int i, int j){
+  FILE *file;
+  file=fopen(ime, "w");
+  if(file==NULL){
+    fprintf(stderr, "Failed opening file. w mode. file=%s\n", ime);
+     exit(1);
+  }
+
+  
+  i=I/2;
+  j=J/2;
+  int k,l;
+  float dir[4];
+  
+  //for(k=0;k<J;k++){
+   //l=i*J*I+j+k*I;
+  for(k=0;k<K;k++){
+   l=i+j*I+k*I*J;
+   calcQ2dirL(l,dir);
+   if((LMARK[l]&LMARKKOL)) fprintf(file,"%d 0 0 0 0\n", k);
+   else fprintf(file,"%d %.5e %.5e %.5e %.5e\n", k,dir[0], dir[1], dir[2], dir[3]);
+   //if(!(LMARK[l]&LMARK1)) fprintf(file,"%d 0 0 0 0\n", k);
+   //else fprintf(file,"%d %.5e %.5e %.5e %.5e\n", k,U[1][l], U[2][l], U[3][l], U[4][l]);
+  }
+  fclose(file);
+}
+
+
+//Writes the velocity field (with Paraview header)
+void zapisUraw(char *ime){
+  // do it in binary
+  int l;
+  float tempU[3];
+  FILE *file;
+
+  file=fopen(ime,"w");
+  if(file==NULL){
+    fprintf(stderr, "Failed opening file. w mode. file=%s\n", ime);
+     exit(1);
+  }
+
+  
+  fprintf(file,"# vtk DataFile Version 3.0\nKnot\nBINARY\nDATASET STRUCTURED_POINTS\n");
+  fprintf(file,"DIMENSIONS %d %d %d\n",I,J,K);
+  fprintf(file,"ORIGIN 0 0 0\n");
+  fprintf(file,"SPACING 1 1 1\n");
+  fprintf(file,"POINT_DATA %d\n",I*J*K);
+  fprintf(file,"VECTORS U float\n");
+
+  for(l=0; l<NMAX; l++)
+  {
+    if(LMARK[l]!=LMARK1)
+    { 
+      tempU[0]=0; tempU[1]=0; tempU[2]=0;
+    }
+    else
+    {
+      // go to bigendian
+      float temp0 = FloatSwap(U[1][l]); 
+      tempU[0]=temp0;
+      float temp1 = FloatSwap(U[2][l]);
+      tempU[1]=temp1;
+      float temp2 = FloatSwap(U[3][l]);
+      tempU[2]=temp2;
+    }
+    fwrite(tempU,sizeof(float),3,file);
+  };
+  fclose(file);
+}
+
+
+// //Writes the density field (with Amira header)
+// void zapisDENraw(char *ime){
+//   
+//   int l;//i,j,k; 
+//   FILE *file;
+//  
+//   file=fopen(ime,"w");
+//   if(file==NULL){
+//     fprintf(stderr, "Failed opening file. w mode. file=%s\n", ime);
+//      exit(1);
+//   }
+// 
+//   fprintf(file,"# AmiraMesh BINARY-LITTLE-ENDIAN 2.1\n\n");
+//   fprintf(file,"define Lattice %d %d %d\n\n",I,J,K);
+//   fprintf(file,"Parameters {\n");
+//   fprintf(file,"BoundingBox 0 %d 0 %d 0 %d\n",I-1,J-1,K-1);
+//   fprintf(file,"CoordType \"uniform\"\n}\n\n");
+//   fprintf(file,"Lattice { float[1] Data } @1\n\n");
+//   fprintf(file,"# Data section follows \n@1\n");
+// 
+//   for(l=0;l<NMAX;l++){ 
+//     float delovni;
+//     /*if(LMARK[l]&LMARKRP){ delovni=DENSITYINIT;}
+//     else*/ delovni=U[4][l];
+//     fwrite(&delovni,sizeof(float),1,file);
+//   }
+//   fclose(file);
+// }
+
+//Writes the density field (with Paraview header)
+void zapisDENraw(char *ime){
+  
+  int l;//i,j,k; 
+  FILE *file;
+  float tempDEN[1];
+ 
+  file=fopen(ime,"w");
+  if(file==NULL){
+    fprintf(stderr, "Failed opening file. w mode. file=%s\n", ime);
+     exit(1);
+  }
+    fprintf(file,"# vtk DataFile Version 3.0\nKnot\nBINARY\nDATASET STRUCTURED_POINTS\n");
+    fprintf(file,"DIMENSIONS %d %d %d\n",I,J,K);
+    fprintf(file,"ORIGIN 0 0 0\n");
+    fprintf(file,"SPACING 1 1 1\n");
+    fprintf(file,"POINT_DATA %d\n",I*J*K);
+    fprintf(file,"SCALARS density float\nLOOKUP_TABLE default\n");
+
+
+  for(l=0;l<NMAX;l++){ 
+    float delovni;
+    /*if(LMARK[l]&LMARKRP){ delovni=DENSITYINIT;}
+    else*/ 
+    float temp4=FloatSwap(U[4][l]);
+    tempDEN[0]=temp4;
+
+    fwrite(tempDEN,sizeof(float),1,file);
+  }
+  fclose(file);
+}
+
+
+
+//Computes the local director field (in the point with index l) from the local Q-tensor
+void calcQ2dirL(int l,float *dir)
+{
+	int i,max_ind;
+      if(!(LMARK[l]&LMARKKOL)){
+	float **del,d[4],e[4],max;
+
+	del=matrix(1,3,1,3);
+
+	//Diagonalisation
+	del[1][1]=Q[1][l]; del[1][2]=Q[4][l]; del[1][3]=Q[5][l];
+	del[2][1]=Q[4][l]; del[2][2]=Q[2][l]; del[2][3]=Q[6][l];
+	del[3][1]=Q[5][l]; del[3][2]=Q[6][l]; del[3][3]=Q[3][l];
+	tred2(del,3,d,e);  //Sintax: matrix to be diagonalised (will be modified in the process); matrix dimension; main diagonal vector; upper diagonal vector
+	tqli(d,e,3,del);   //Sintax: main diagonal (where eigenvalues are later stored); upper+bottom diagonal; dimension; output matrix from tred2
+
+	//Find the largest eigenvalue and the director
+	max=d[1]; max_ind=1;
+	for(i=2;i<=3;i++){if(d[i]>max){max=d[i];max_ind=i;};};
+	dir[0]=del[1][max_ind];dir[1]=del[2][max_ind];dir[2]=del[3][max_ind];
+	dir[3]=max;
+	//printf("%f %f %f     ", dir[1], dir[2], dir[3]);
+	free_matrix(del,1,3,1,3);
+      }
+      else{ for(i=0;i<3;i++) dir[i]=0.f; dir[3]=1;}
+}
+
+
+//Writes the director field (with Amira header), computed from the Q-tensor
+void zapisQ2DIRraw(char *ime)
+{
+	int l;
+	float dir[4];
+	FILE *file;
+
+// 	file=fopen(ime,"w");
+// 	fprintf(file,"# AmiraMesh BINARY-LITTLE-ENDIAN 2.1\n\n");
+// 	fprintf(file,"define Lattice %d %d %d\n\n",I,J,K);
+// 	fprintf(file,"Parameters {\n");
+// 	fprintf(file,"BoundingBox 0 %d 0 %d 0 %d\n",I-1,J-1,K-1);
+// 	fprintf(file,"CoordType \"uniform\"\n}\n\n");
+// 	fprintf(file,"Lattice { float[3] Data } @1\n\n");
+// 	fprintf(file,"# Data section follows \n@1\n");
+
+    
+    file=fopen(ime,"w");
+    if(file==NULL){
+        fprintf(stderr, "Failed opening file. w mode. file=%s\n", ime);
+         exit(1);
+    }
+
+    fprintf(file,"# vtk DataFile Version 3.0\nKnot\nBINARY\nDATASET STRUCTURED_POINTS\n");
+    fprintf(file,"DIMENSIONS %d %d %d\n",I,J,K);
+    fprintf(file,"ORIGIN 0 0 0\n");
+    fprintf(file,"SPACING 1 1 1\n");
+    fprintf(file,"POINT_DATA %d\n",I*J*K);
+    fprintf(file,"VECTORS n float\n");
+
+	for(l=0; l<NMAX; l++){
+	  calcQ2dirL(l,dir);														//calculate director locally
+	  //if(LMARK[l]&LMARK1){dir[0]=0; dir[1]=0; dir[2]=0;}
+	  //else{dir[0]=QNIC[1][l]; dir[1]=QNIC[2][l]; dir[2]=QNIC[3][l];}
+     
+      
+      // go to bigendian
+      float dir0 = FloatSwap(dir[0]); 
+      dir[0]=dir0;
+      float dir1 = FloatSwap(dir[1]);
+      dir[1]=dir1;
+      float dir2 = FloatSwap(dir[2]);
+      dir[2]=dir2;
+      
+      fwrite(dir,sizeof(float),3,file);
+        
+    };
+    fclose(file);
+        
+    
+}
+
+
+//Writes the order parameter field (with Amira header), computed from the Q-tensor
+void zapisQ2OPraw(char *ime)
+{
+	/*int l;
+	float dir[4];
+	float delovni;
+	FILE *file;
+
+	file=fopen(ime,"w");
+	fprintf(file,"# AmiraMesh BINARY-LITTLE-ENDIAN 2.1\n\n");
+	fprintf(file,"define Lattice %d %d %d\n\n",I,J,K);
+	fprintf(file,"Parameters {\n");
+	fprintf(file,"BoundingBox 0 %d 0 %d 0 %d\n",I-1,J-1,K-1);
+	fprintf(file,"CoordType \"uniform\"\n}\n\n");
+	fprintf(file,"Lattice { float[1] Data } @1\n\n");
+	fprintf(file,"# Data section follows \n@1\n");
+
+	for(l=0; l<NMAX; l++){
+	  calcQ2dirL(l,dir);														//calculate director locally
+	  delovni=dir[3];
+	  fwrite(&delovni,sizeof(float),1,file);
+	};
+	fclose(file);*/
+
+    
+    
+    //Writes the order parameter field (with Paraview header), computed from the Q-tensor    
+    // do it in binary
+    int l;
+    float dir[4];
+    float delovni;
+    FILE *file;
+    
+    float Sk=0.4;
+    float S;
+    float DFV=0;
+
+    file=fopen(ime,"w");
+    if(file==NULL){
+        fprintf(stderr, "Failed opening file. w mode. file=%s\n", ime);
+        exit(1);
+    }
+
+    fprintf(file,"# vtk DataFile Version 3.0\nKnot\nBINARY\nDATASET STRUCTURED_POINTS\n");
+    fprintf(file,"DIMENSIONS %d %d %d\n",I,J,K);
+    fprintf(file,"ORIGIN 0 0 0\n");
+    fprintf(file,"SPACING 1 1 1\n");
+    fprintf(file,"POINT_DATA %d\n",I*J*K);
+    fprintf(file,"SCALARS order float\nLOOKUP_TABLE default\n");
+    
+//     FILE *defect_volume;
+//     defect_volume=fopen("./data/defect_volume.txt","a");
+//     fprintf(defect_volume, "DFV\n");
+
+    int count=0;
+    for(l=0; l<NMAX; l++){
+        calcQ2dirL(l,dir);														//calculate director locally    
+        delovni = FloatSwap(dir[3]); // go to bigendian
+        fwrite(&delovni,sizeof(float),1,file);
+        
+        S=dir[3];
+        if (LMARK[l]==LMARK1) {
+            count++;
+            if (S < Sk) DFV+=1;
+        }
+
+        
+         
+        
+    };
+    DFV=DFV/count;
+//     fprintf(defect_volume, "%f\n", DFV);
+//     fclose(defect_volume);
+//     printf("defect volume fraction: %f\n", DFV);
+    
+    fclose(file);
+//     fclose(defect_volume);
+    
+
+}
+
+
+//Sum of points where S is smaller than Sk=0.4
+
+void zapisDVF(char *ime){
+    int l;
+    float dir[4];
+    float delovni;
+    
+    float Sk=0.4;
+    float S;
+    float DVF=0;
+    
+    FILE *defect_volume;
+    defect_volume=fopen("./data/defect_volume.txt","a");
+    if(defect_volume==NULL){
+        fprintf(stderr, "Failed opening file. a mode. file=%s\n", "./data/defect_volume.txt");
+        exit(1);
+    }
+
+
+
+    int count=0;
+    for(l=0; l<NMAX; l++){
+        calcQ2dirL(l,dir);														//calculate director locally    
+        delovni = FloatSwap(dir[3]); // go to bigendian
+        S=dir[3];
+        if (LMARK[l]==LMARK1) {
+            count++;
+            if (S < Sk) DVF+=1;
+        }
+        
+    };
+    
+    DVF=DVF/count;
+    fprintf(defect_volume, "%f\n", DVF);
+//     fclose(defect_volume);
+    printf("defect volume fraction: %f\n", DVF);
+    fprintf(stderr,"defect volume fraction: %f\n", DVF);
+
+    fclose(defect_volume);
+    
+}
+
+//----------------------------------------Dipol distance from coloid------------------------------------------------------------------------
+
+//Points where S is smaller than Sk=0.4 and theirs distance from coloid's center (0,0,0)
+void zapisDistance(){
+    float dir[4];
+    float delovni;
+    
+    FILE *dipol_distance;
+    dipol_distance=fopen("./data/dipol_distance.txt","a");
+    if(dipol_distance==NULL){
+        fprintf(stderr, "Failed opening file. a mode. file=%s\n", "./data/dipol_distance.txt");
+        exit(1);
+    }
+
+        
+    float Sk=0.4;
+    float S;
+    float DVF=0;
+    float distance=0;
+    float second_distance=0;    
+    int count=0;
+    double center[3]={0,0,0};
+//     fprintf(stderr, "sem v funkciji dipol distance\n");
+     
+    
+    for(int l=0; l<NMAX; l++){
+        calcQ2dirL(l,dir);														//calculate director locally    
+        delovni = FloatSwap(dir[3]); // go to bigendian
+//         printf("for zanka v dipol distance\n");
+        S=dir[3];
+
+//         fprintf(stderr, "LMARK[l] = %d; ",LMARK[l]);
+//         
+        if(LMARK[l]==LMARK1){
+//             fprintf(stderr, "S = %f; ",S);
+            if (S < Sk){
+                count++;
+                float x,y,z,r,rr; 
+                            
+                x=(float)(i_vr(l))-(float)(I-1)/2+0.001;
+                y=(float)(j_vr(l))-(float)(J-1)/2+0.001;
+                z=(float)(k_vr(l))-(float)(K-1)/2+0.001;
+                
+                rr = (x - center[0]) * (x - center[0]) + (y - center[1]) * (y - center[1]) + (z - center[2]) * (z - center[2]);
+                r = sqrt((x - center[0]) * (x - center[0]) + (y - center[1]) * (y - center[1]) + (z - center[2]) * (z - center[2]));
+                distance = distance + r;
+//                 fprintf(stderr, "distance = %f; ",distance);
+//                 fprintf(stderr, "count = %d\n",count);
+                second_distance = second_distance + rr;
+                
+            
+            }
+            
+        };
+        
+
+
+    };
+    
+//     fprintf(stderr, "konec funkcije\n");
+  
+    fprintf(stderr, "%f,%f,%f,%f\n",distance, distance/count, second_distance, second_distance/count);
+    fprintf(dipol_distance, "%f,%f,%f,%f\n",distance, distance/count, second_distance, second_distance/count);
+    fclose(dipol_distance);
+    
+}
+
+
+//-----------------------------------------------DIMERI---------------------------------------------------------------------------------------
+// //Reads the data file for active droplet
+// void beriDROPLETDATA(char *ime)
+// {
+// 	FILE *file;
+// 	int l,m;
+// 	float delovni;
+// 	char delovni2;
+// 	int result;
+// 
+// 	file=fopen(ime,"r");
+//     if(file==NULL){
+//         fprintf(stderr, "Failed opening file. r mode. file=%s\n", ime);
+//          exit(1);
+//     }
+// 
+//     fprintf(stderr,"zacetek\n");
+//     
+//     int sizes[3]; 
+//     int q0;
+//     int count=0;
+//     float normal[3];
+//     fread(sizes,sizeof(int),3,file);
+//     fprintf(stderr,"beremo velikost %dx%dx%d, Nmax=%d\n",sizes[0],sizes[1],sizes[2],NMAX);
+//     
+//       
+//     
+// 	for(l=0;l<NMAX;l++){
+// 		for(m=0;m<6;m++)
+// 		{
+// 		  result=fread(&delovni,sizeof(float),1,file);
+// 		  if (result != 1) {fputs ("Error reading data!\n",stderr); exit (2);}
+// 		  Q[m+1][l]=delovni;
+//           QNEW[m+1][l]=delovni;
+// 		};
+//         
+//         UNIC[1][l]=0.f;		//surface velocity
+//         UNIC[2][l]=0.f;
+//         UNIC[3][l]=0.f;
+//         
+//         
+//         //change LMARK      
+//         result=fread(&delovni2,sizeof(char),1,file);
+//         
+//         LMARK[l]=delovni2;
+//         if (delovni2==1) {LMARK[l]=32;}
+//         else if (delovni2==16 || delovni2==32 || delovni2==64 || delovni2==96) {LMARK[l]=8;}
+// /*        delovni2=LMARK[l];
+//         fwrite(&delovni2,sizeof(char),1,file); */        
+//         
+//         
+// //         if(delovni2)fprintf(stderr,"%d %d\n",l,delovni2);
+//         //if(l%3==0)        fprintf(stderr,"%f %f %f %f %f %f\n",Q[1][l],Q[2][l],Q[3][l],Q[4][l],Q[5][l],Q[6][l]);
+// 	}
+// 	    for(l=0;l<NMAX;l++){
+//         int lxp,lxm,lyp,lym,lzp,lzm;
+//         if(i_vr(l)==0) {lxp=l+1; lxm=l+I-1;}		//periodicnost po x smeri
+//         else if(i_vr(l)==I-1) {lxp=l-(I-1); lxm=l-1;}
+//         else {lxp=l+1; lxm=l-1;};
+// 
+//         if(j_vr(l)==0) {lyp=l+I; lym=l+(J-1)*I;}	//periodicnost po y smeri
+//         else if(j_vr(l)==J-1) {lyp=l-(J-1)*I; lym=l-I;}
+//         else {lyp=l+I; lym=l-I;};
+// 
+//         if(k_vr(l)==0) {lzp=l+I*J; lzm=l+(K-1)*I*J;}	//periodicnost po z smeri
+//         else if(k_vr(l)==K-1) {lzp=l-(K-1)*I*J; lzm=l-I*J;}
+//         else {lzp=l+I*J; lzm=l-I*J;}; 
+//         
+//         
+//         if(LMARK[lxp]&(LMARKPOV|LMARKKOL) && LMARK[lxm]&(LMARKPOV|LMARKKOL) && LMARK[lyp]&(LMARKPOV|LMARKKOL) && LMARK[lym]&(LMARKPOV|LMARKKOL) && LMARK[lzp]&(LMARKPOV|LMARKKOL) && LMARK[lzm]&(LMARKPOV|LMARKKOL)){ LMARK[l]=(LMARKKOL|LMARK0);}
+//   }
+// 
+//   for(l=0;l<NMAX;l++){
+//     if(i_vr(l)==0 || i_vr(l)==I-1 || j_vr(l)==0 || j_vr(l)==J-1 || k_vr(l)==0 || k_vr(l)==K-1) LMARK[l]=(LMARK[l]|LMARKRP);														
+//     else if(LMARK[l]&LMARKKOL) for(int m=1; m<ST_HITROSTI; m++) if(LMARK[calcLBlnew(l, m)]&(LMARK1|LMARKPOV)) LMARK[l]=LMARKKOL;
+//   }
+//   
+// 	
+// 	result=fread(&q0, sizeof(int),1,file);
+//     fprintf(stderr,"result = %d\n",result);
+//     
+//     fprintf(stderr,"q0 = %d\n",q0);
+//     for(l=0;l<NMAX;l++){
+//         if(LMARK[l]==8){
+//             fread(normal, sizeof(float), 3, file);
+//             count++;
+//             QNIC[1][l]= normal[0];
+//             QNIC[2][l]= normal[1];
+//             QNIC[3][l]= normal[2];
+//             /*
+//             
+//             QNIC[1][l]=0.5f*SEQ*(3*normal[0]*normal[0]-1);
+//             QNIC[2][l]=0.5f*SEQ*(3*normal[1]*normal[1]-1);
+//             QNIC[3][l]=0.5f*SEQ*(3*normal[2]*normal[2]-1);
+//             QNIC[4][l]=1.5f*SEQ*(normal[0]*normal[1]);
+//             QNIC[5][l]=1.5f*SEQ*(normal[0]*normal[2]);
+//             QNIC[6][l]=1.5f*SEQ*(normal[1]*normal[2]);
+//             */
+//         }
+//     }
+//     
+//     fprintf(stderr,"count = %d\n",count);
+// 	
+//     //fprintf(stderr,"konec\n");
+// 	fclose(file);
+// }  
+
+
+

dipol/fibonacci/saturn_ring/state_A1_test/acn_utils.c

@@ -0,0 +1,81 @@
+
+
+//Line number in the x direction. Range is from 0 to I-1.
+ int i_vr(int l){
+ return l%I; 
+}
+
+//Line number in the y direction
+ int j_vr(int l){
+ return (l%(I*J))/I; 
+}
+
+//Line number in the z direction
+ int k_vr(int l){
+ return l/(I*J); 
+}
+
+//Characteristic velocity vectors of the LV model. Opposite vector pairs have to be written in odd-even combination (important for calcLBmnew function)!
+void initialiseE()
+{
+  E[0][1]= 0;		E[0][2]= 0;		E[0][3]= 0;
+  
+  E[1][1]= 1;		E[1][2]= 0;		E[1][3]= 0;
+  E[2][1]= -1;		E[2][2]= 0;		E[2][3]= 0;
+  E[3][1]= 0;		E[3][2]= 1;		E[3][3]= 0;
+  E[4][1]= 0; 		E[4][2]= -1;		E[4][3]= 0;
+  E[5][1]= 0;		E[5][2]= 0;		E[5][3]= 1;
+  E[6][1]= 0;  		E[6][2]= 0;  		E[6][3]= -1;
+  
+  E[7][1]= 1;  		E[7][2]= 1;	  	E[7][3]= 0;
+  E[8][1]= -1; 		E[8][2]= -1;		E[8][3]= 0;
+  E[9][1]= -1;  	E[9][2]= 1;		E[9][3]= 0;
+  E[10][1]= 1;		E[10][2]= -1;		E[10][3]= 0;
+  
+  E[11][1]= 1;  	E[11][2]= 0;		E[11][3]= 1;
+  E[12][1]= -1;  	E[12][2]= 0;  		E[12][3]= -1;
+  E[13][1]= -1;  	E[13][2]= 0;		E[13][3]= 1;
+  E[14][1]= 1;		E[14][2]= 0;		E[14][3]= -1;
+  
+  E[15][1]= 0;  	E[15][2]= 1;		E[15][3]= 1;
+  E[16][1]= 0;  	E[16][2]=-1;  		E[16][3]= -1;
+  E[17][1]= 0;  	E[17][2]= -1;		E[17][3]= 1;
+  E[18][1]= 0;		E[18][2]= 1;		E[18][3]= -1;
+  
+  //Constants for the model. Adapted from Thurey, Rude, Korner
+  int i;
+  WKONST[0]=1.f/3.f;
+  for(i=1;i<=6; i++) WKONST[i]=1.f/18.f;
+  for(i=7;i<=18; i++) WKONST[i]=1.f/36.f;
+}
+
+
+
+//Function that multiplies matrices Aik x Bkj = Cij t and returns the element Cij
+void AMatrixB(float (*c)[4], float (*a)[4],float (*b)[4])
+{
+  for (int i = 1; i<=3; i++)
+    {
+      for (int j = 1; j<=3; j++)
+	{
+	  c[i][j]=a[i][1]*b[1][j]+a[i][2]*b[2][j]+a[i][3]*b[3][j];
+	}
+    }
+}
+
+// Function to take a little-endian encoded float and make it big-endian (or vice versa). Useful for reading to paraview
+float FloatSwap( float f )
+{
+    union
+    {
+        float f;
+        char b[4];
+    } dat1, dat2;
+
+    dat1.f = f;
+    dat2.b[0] = dat1.b[3];
+    dat2.b[1] = dat1.b[2];
+    dat2.b[2] = dat1.b[1];
+    dat2.b[3] = dat1.b[0];
+    return dat2.f;
+}

dipol/fibonacci/saturn_ring/state_A1_test/an_01.c

@@ -0,0 +1,218 @@
+//This is a hybrid lattice-Boltzmann code for nematic hydrodynamics
+
+//Create folder data!!
+//Compile with "gcc -O3 -o prog an_01.c -lm -lpthread -march=native -std=gnu99 -fopenmp -Wall -Wno-unused-variable -Wno-unused-but-set-variable"
+//Compile with "gcc -O3 -o prog an_01.c -lm -lpthread -march=native -std=gnu99 -fopenmp -Wall -Wno-unused-variable -Wno-unused-but-set-variable -Wno-format-overflow"
+
+//ActiveDroplet 
+//Compile with "gcc -O3 -o prog2 an_01.c diff_eq_system.c -lm -lpthread -march=native -std=gnu99 -fopenmp -Wall -Wno-unused-variable -Wno-unused-but-set-variable -Wno-format-overflow"
+//gcc -O3 -o prog2 -g an_01.c diff_eq_system.c -lm -lpthread -march=native -std=gnu99 -fopenmp -Wall -Wno-unused-variable -Wno-unused-but-set-variable -Wno-format-overflow
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <omp.h>
+#include <string.h>
+
+
+int I;	//Number of points in each direction
+int J;
+int K;
+#define NMAX			(I*J*K)		//NMAX is no longer defined as (I+1)*(J+1)*(K+1)
+
+#define ST_HITROSTI		19		//Number of velocity vectors within the model -- in this case 3DQ19 -- 19 velocities		To change the model, one needs to set new values of ST_HITROSTI, new vectors E and also set FEQ, P, and bounce-back
+
+#define STCASKOR		5000000
+#define STCASKORINIT		0
+
+#define STCASKORQ		1		//How many steps in Q are performed for each step in LB
+
+#define STKONV			100		//On how many steps the fields are written to a file
+#define STKONVV			10000		//On how many steps the fields are written to a file
+#define STKONG          100     //On how many steps the defect grad and config are written to a file
+
+
+#define TWRITEU			1		//Writes velocity
+#define TWRITEDEN		1		//Writes density
+#define TWRITEDIR		1		//Writes director field
+#define TWRITEOP		1       //Writes order parameter field
+#define TWRITEQ			1       //Writes Q-tensor field
+
+#define STDVF			100		//On how many steps the defect volume fraction is written to a file
+#define TWRITEDVF		1       //Writes defect volume fraction
+#define STFE			500		//On how many steps the free energy is written to a file
+#define TWRITEFE		1       //Writes free energy
+#define TWRITELOC       1       //Writes location of defects
+#define TWRITEgrad      1       //Writes defet gradients
+
+// #define STU             300     //Change the initial condition
+// #define initialU        1
+
+#define STCOPY          10000      //On how many steps the copy of DATA.bin is written to a file 
+#define TWRITECOPY      1      //Writes copy of DATA.bin
+
+
+
+#define ZAPISF			1		//Write all necessary functions to a file at the end of the computation
+#define BERIF			0		//Read initialisation from the file
+
+#define READINITIAL     0
+
+
+
+
+// #define ZETAchange		70	    //change zeta value
+// 
+// #define ZETAstart		50	//change zeta value (start)
+// #define ZETAend		    80	//change zeta value (end)
+// 
+
+int STPROC=64;					//Number of CPU cores for the paralelisation
+
+float **U, **F, **FEQ, **FNEW, **P;				//Allocation of global matrices
+float **Q, **QNIC, **QNEW;
+float *WKONST;
+int **E;
+char *LMARK;
+float **UNIC;
+float **H;
+float **FORCE;
+float **TAU;
+
+char* QMARK;
+
+// number of surface data points
+int Qmax;
+float* QNC;
+
+float **gradU;
+float **gradU_defect;
+
+char Izpis[] = "droplet";
+char imeIzpis[100];
+
+
+#include "acn_funk_head.h"		//function headers
+#include "diff_eq_system.h"
+#include "konstante.c"
+#include "num_rec.c"			//All needed function from numerical recipes in a single file
+#include "acn_utils.c"			//Useful functions -- utility
+#include "acn_active_droplet.c"    //defects on shell
+#include "acn_funk_lb.c"		//Lattice-Boltzmann part of the code
+#include"acn_funk_op.c"
+#include "acn_funk_zapis.c"
+
+
+
+
+int main(int argc, char** args){
+
+  ZETA=(float) atof(args[1]);
+  I=(int) atoi(args[2]);
+  J=(int) atoi(args[3]);
+  K=(int) atoi(args[4]);
+  
+  sprintf(imeIzpis, "%s-%dx%dx%d_%f", Izpis, I, J, K, ZETA);
+  
+  long Seed=12345; 					//Seed for random initial conditions
+  char buf[256];
+  int l,m;
+  
+  U=matrix(1,4,0,NMAX-1);				// velocity field; 4th component is density
+  F=matrix(0,ST_HITROSTI-1,0,NMAX-1);		// Lattice Boltzmann scalar distribution functions to calculate velocity U
+  FNEW=matrix(0,ST_HITROSTI-1,0,NMAX-1);		// matrix field to temporarily store currently calculated F
+  FEQ=matrix(0,ST_HITROSTI-1,0,NMAX-1);		// Lattice Boltzmann scalar distribution functions to calculate velocity U - equilibrium; nedded in colloision operators
+  P=matrix(0,ST_HITROSTI-1,0,NMAX-1);		// Lattice Boltzmann scalar distribution functions to calculate velocity U - forcing term; needed in colloision operators
+  TAU=matrix(1,9,0,NMAX-1);
+  
+  Q=matrix(1,6,0,NMAX-1);				//Q-tensor
+  QNEW=matrix(1,6,0,NMAX-1);
+  QNIC=matrix(1,6,0,NMAX-1);			//Qnic include info about the normal, prefered surface orientation and possibly also surface velocity
+  
+  WKONST=vector(0,18);
+  LMARK=malloc(NMAX*sizeof(char));
+  UNIC=matrix(1,3,0,NMAX-1);			//Surface velocity (usually commented out to save time and RAM)
+  gradU=matrix(1,3,0,NMAX-1);
+  gradU_defect=matrix(1,3,0,NMAX-1);
+  
+  H=matrix(1,6,0,NMAX-1);				//H tensor -- molecular field
+  FORCE=matrix(1,3,0,NMAX-1);				//FORCE = divergence of a total stress tensor -- needed to correct the velocity due to discrete effects
+  
+  E=imatrix(0,ST_HITROSTI-1,1,3);			// matrix of all LB characteristic velocity lattice vectors
+  
+  initialiseE();				//Characteristic lattice vectors for the chosen LB model
+
+// if(BERIF){ sprintf(buf, "droplet-200x200x200_0.010000-DATA.bin", imeIzpis); beriDATA(buf);		//Reads F
+  
+  if(BERIF){ sprintf(buf, "%s-DATA.bin", imeIzpis); beriDATA(buf);		//Reads F
+    for(l=0;l<NMAX;l++){								
+      for(m=0;m<ST_HITROSTI;m++) FNEW[m][l]=F[m][l];				//Transcribe F to FNEW
+      FORCE[1][l]=0; FORCE[2][l]=0; FORCE[3][l]=0; calcF2Ul(l,F);		//Compute U
+    }
+  }
+  else{
+    zacPriblQ(Seed);
+    zacPriblUF(Seed);				// initial condition for U and F
+   
+  }
+  
+  //read initial conditions (janus core)
+//   if(READINITIAL) {sprintf(buf,"dim_par_eightL_relax_E01.raw"); beriDROPLETDATA(buf);} 
+//   if(READINITIAL) {sprintf(buf,"dim_par_satrn_relax_E01.raw"); beriDROPLETDATA(buf);}
+
+
+//     ActiveDroplet sphere1;
+//     ActiveDroplet* sphere1_ptr = &sphere1;
+//     init_droplet(sphere1_ptr);
+//     printf("Len after init: %d\n", sphere1.num_points);
+//printf("%d\n", sphere1.num_points);
+      
+  izrQU();					//Main computational loop
+
+  
+  
+  //Calculates the upper estimates for the Reynolds and Ericksen number
+  float umax=0;
+  for(l=0;l<NMAX;l++) if(LMARK[l]&LMARK1){
+   float delovni=fsqr(U[1][l]) + fsqr(U[2][l]) + fsqr(U[3][l]);
+   if(delovni>umax) umax=delovni;
+  }
+  printf("reynolds: %e\nericksen: %e\n gostota: %f\n", 3*sqrt(umax)*I/(TAUF/DT-0.5f), 1.f/3.f*(TAUF/DT-0.5f)*DENSITYINIT*sqrt(umax)*I, DENSITYINIT);
+  
+
+  
+  
+  //Write the desired variables in .raw files
+  //Write the desired variables in .vtk files
+  sprintf(buf,"./data/%s_u.vtk",imeIzpis);
+  zapisUraw(buf);					// write velocity profile
+
+  sprintf(buf,"./data/%s_den.vtk",imeIzpis);
+  zapisDENraw(buf);					// write density profile
+  
+  sprintf(buf,"./data/%s_dir.vtk",imeIzpis);
+  zapisQ2DIRraw(buf);
+  
+  sprintf(buf,"./data/%s_op.vtk",imeIzpis);
+  zapisQ2OPraw(buf);
+  
+  //sprintf(buf,"./data/%s_u.dat",imeIzpis);		//velocity profile for gnuplot
+  //za_profilU(buf,K/2,I/2);
+  
+//   zapisDVF(buf);
+  
+
+
+  //Write all the variables in a .bin file
+  if(ZAPISF){ sprintf(buf, "%s-DATA.bin", imeIzpis); zapisDATA(buf);}
+  
+  
+  //free all matrices	
+  free_matrix(U,1,4,0,NMAX-1);			free_matrix(F,0,ST_HITROSTI-1,0,NMAX-1);		free_matrix(FNEW,0,ST_HITROSTI-1,0,NMAX-1);	free_matrix(FEQ,0,ST_HITROSTI-1,0,NMAX-1);			free_matrix(P,0,ST_HITROSTI-1,0,NMAX-1);	free_imatrix(E,0,ST_HITROSTI-1,1,3);			free_matrix(Q,1,6,0,NMAX-1);	free_matrix(QNEW,1,6,0,NMAX-1);
+  free_matrix(QNIC,1,6,0,NMAX-1);		free_vector(WKONST,0,18);				free(LMARK);					free_matrix(UNIC,1,3,0,NMAX-1);
+  free_matrix(TAU,1,9,0,NMAX-1);
+  free_matrix(H,1,6,0,NMAX-1);			free_matrix(FORCE,1,3,0,NMAX-1); free_matrix(gradU,1,3,0,NMAX-1); free_matrix(gradU_defect,1,3,0,NMAX-1);
+
+  return 0;
+}

dipol/fibonacci/saturn_ring/state_A1_test/diff_eq_system.c

@@ -0,0 +1,1097 @@
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <complex.h>
+#include <omp.h>
+#include "diff_eq_system.h"
+
+
+
+void ephi(double phi, double* vector) {
+  vector[0] = -sin(phi);
+  vector[1] = cos(phi);
+  vector[2] = 0;
+}
+
+
+void etheta(double phi, double theta, double* vector) {
+  vector[0] = cos(phi) * cos(theta);
+  vector[1] = sin(phi) * cos(theta);
+  vector[2] = -sin(theta);
+}
+
+void sph_to_cart(double theta, double phi, double* cartesian) {
+  cartesian[0] = cos(phi) * sin(theta);
+  cartesian[1] = sin(phi) * sin(theta);
+  cartesian[2] = cos(theta);
+}
+
+
+void cross_product(double* vector0, double* vector1, double* result) {
+
+  result[0] = vector0[1] * vector1[2] - vector0[2] * vector1[1];
+  result[1] = vector0[2] * vector1[0] - vector0[0] * vector1[2];
+  result[2] = vector0[0] * vector1[1] - vector0[1] * vector1[0];
+}
+
+
+void lattice_base(int n, double* thetas, double* phis, double* basis) {
+    // Returns n x 3 x 3 array of local lattice bases, serves as transformation from global to local basis
+    double cos_phi, sin_phi, cos_theta, sin_theta;
+    for (int i=0; i<n; i++) {
+      cos_phi = cos(phis[i]);
+      sin_phi = sin(phis[i]);
+      cos_theta = cos(thetas[i]);
+      sin_theta = sin(thetas[i]);
+      basis[9*i + 0] = cos_phi * sin_theta;
+      basis[9*i + 1] = sin_theta * sin_phi;
+      basis[9*i + 2] = cos_theta;
+      basis[9*i + 3] = cos_theta * cos_phi;
+      basis[9*i + 4] = cos_theta * sin_phi;
+      basis[9*i + 5] = -sin_theta;
+      basis[9*i + 6] = -sin_phi;
+      basis[9*i + 7] = cos_phi;
+      basis[9*i + 8] = 0;
+    }
+}
+
+void rodrigues_rotation(double* vec, double angle, double* axis, double* result) {
+  double cos_angle = cos(angle);
+  double sin_angle = sin(angle);
+  double scalar;
+  double cross[3];
+
+  cross_product(axis, vec, cross);
+  scalar = axis[0]*vec[0] + axis[1]*vec[1] + axis[2]*vec[2];
+
+  result[0] = vec[0] * cos_angle + cross[0] * sin_angle + axis[0] * scalar * (1 - cos_angle);
+  result[1] = vec[1] * cos_angle + cross[1] * sin_angle + axis[1] * scalar * (1 - cos_angle);
+  result[2] = vec[2] * cos_angle + cross[2] * sin_angle + axis[2] * scalar * (1 - cos_angle);
+}
+
+
+void x_to_z_transf(int n, double* vectors) {
+  // rotation of all vectors in array for pi/2 around y-axis (transformation x -> z and z -> -x)
+  double new_x;
+  double new_z;
+
+  for (int i=0; i<n;   i++) {
+    new_x = -vectors[3*i+2];
+    new_z = vectors[3*i];
+    vectors[3*i] = new_x;
+    vectors[3*i+2] = new_z;
+  }
+
+}
+
+
+void inverse_x_to_z_transformation(int n, double* vectors) {
+  // rotation of all vectors in array for -pi/2 around y-axis (transformation x -> -z and z -> x)
+  double new_x;
+  double new_z;
+
+  for (int i=0; i<n;   i++) {
+    new_x = vectors[3*i+2];
+    new_z = -vectors[3*i];
+    vectors[3*i] = new_x;
+    vectors[3*i+2] = new_z;
+  }
+
+}
+
+
+double complex stereogr_proj(double theta, double phi) {
+  return cexp(I * phi) / tan(theta / 2.);
+}
+
+
+void psi_avg(double* thetas, double* phis, double alpha0, double* result) {
+  // Minimizing defect orientation
+  double complex ln_sum;
+  for (int i=0; i<4; i++) {
+    ln_sum = 0;
+    for (int j=0; j<4; j++) {
+      if (j != i) {
+        ln_sum += clog(stereogr_proj(thetas[i], phis[i]) - stereogr_proj(thetas[j], phis[j]));
+      }
+    }
+    result[i] = phis[i] - 2 * alpha0 - cimag(ln_sum);
+  }
+}
+
+
+double cos_beta_ij(double theta_i, double phi_i, double theta_j, double phi_j) {
+  // Angular distance between defects
+  return cos(theta_i) * cos(theta_j) + sin(theta_i) * sin(theta_j) * cos(phi_i - phi_j);
+}
+
+
+double dcos_beta_dtheta(double theta_i, double phi_i, double theta_j, double phi_j) {
+  // d(cos beta_ki) / d theta_k
+  return -sin(theta_i) * cos(theta_j) + cos(theta_i) * sin(theta_j) * cos(phi_i - phi_j);
+}
+
+
+double dcos_beta_dphi(double theta_i, double phi_i, double theta_j, double phi_j) {
+  // d(cos beta_ki) / d phi_k
+  return  -sin(theta_i) * sin(theta_j) * sin(phi_i - phi_j);
+}
+
+
+double free_energy(double* thetas, double* phis, double* delta_psi, double chi) {
+  // elastic free energy
+  double en = 0;
+  for (int i=0; i<4; i++) {
+    en += chi * delta_psi[i] * delta_psi[i];
+    for (int j=0; j<4; j++) {
+      if (j != i) {
+        en += -clog(1 - cos_beta_ij(thetas[i], phis[i], thetas[j], phis[j]));
+      }
+    }
+  }
+  return en;
+}
+
+
+double free_energy_pos(double* thetas, double* phis) {
+  // lastic free energy
+  double en = 0;
+  for (int i=0; i<4; i++) {
+    for (int j=0; j<4; j++) {
+      if (j != i) {
+        en += -clog(1 - cos_beta_ij(thetas[i], phis[i], thetas[j], phis[j]));
+      }
+    }
+  }
+  return en;
+}
+
+
+double free_energy_orient(double* delta_psi, double chi) {
+  // lastic free energy
+  double en = 0;
+  for (int i=0; i<4; i++) {
+    en += chi * delta_psi[i] * delta_psi[i];
+  }
+  return en;
+}
+
+
+void df_dtheta(double* thetas, double* phis, double* delta_psi, double chi, double* derivatives) {
+  // Free energy derivative over theta_k
+  for (int k=0; k<4; k++) {
+    derivatives[k] = 0;
+    for (int i=0; i<4; i++) {
+      if (i != k) {
+        derivatives[k] += 2 * (
+                          dcos_beta_dtheta(thetas[k], phis[k], thetas[i], phis[i]) /
+                          (1 - cos_beta_ij(thetas[k], phis[k], thetas[i], phis[i])) -
+                          chi * (delta_psi[k] + delta_psi[i]) / sin(thetas[k]) *
+                          cimag(stereogr_proj(thetas[k], phis[k]) /
+                          (stereogr_proj(thetas[k], phis[k]) - stereogr_proj(thetas[i], phis[i])))
+                        );
+      }
+    }
+  }
+}
+
+
+void df_dphi(double* thetas, double* phis, double* delta_psi, double chi, double* derivatives) {
+  // Free energy derivative over phi_k
+  for (int k=0; k<4; k++) {
+    derivatives[k] = 0;
+    for (int i=0; i<4; i++) {
+      if (i != k) {
+        derivatives[k] += 2 * (
+                          dcos_beta_dphi(thetas[k], phis[k], thetas[i], phis[i]) /
+                          (1 - cos_beta_ij(thetas[k], phis[k], thetas[i], phis[i])) +
+                          chi * (
+                            (delta_psi[k] + delta_psi[i]) *
+                            creal(stereogr_proj(thetas[k], phis[k]) /
+                            (stereogr_proj(thetas[k], phis[k]) - stereogr_proj(thetas[i], phis[i]))) -
+                            delta_psi[k] / 3.
+                          )
+                        );
+      }
+    }
+  }
+}
+
+
+double df_dalpha0(double* delta_psi, double chi) {
+  // Free energy derivative over alpha0
+  double derivative = 0;
+  for (int k=0; k<4; k++) {
+    derivative += 4 * chi * delta_psi[k];
+  }
+  return derivative;
+}
+
+
+void df_dpsi(double* delta_psi, double chi, double* derivatives) {
+  // Free energy derivative over psi_k
+  for (int k=0; k<4; k++) {
+    derivatives[k] = 2 * chi * delta_psi[k];
+  }
+}
+
+
+void eq_system(double* variables, double nu, double chi, double xi_rot, double xi_alpha, double* dy_dt) {
+    // System of differential equations. y=[theta_k, phi_k, psi_k, alpha0] (len(y)=13)
+    double thetas[4];
+    double phis[4];
+    double psis[4];
+    double delta_psi[4];
+    double alpha0 = variables[12];
+    double average_psi[4];
+
+    psi_avg(&variables[0], &variables[4], alpha0, average_psi);
+
+    for (int i=0; i<4; i++) {
+      thetas[i] = variables[i];
+      phis[i] = variables[4 + i];
+      psis[i] = variables[8 + i];
+      delta_psi[i] = psis[i] - average_psi[i];
+      delta_psi[i] -= M_PI * (2 * floor((delta_psi[i]+M_PI) / (2 * M_PI)));  // map into the interval [-pi , pi]
+    }
+
+    double theta_deriv[4];
+    double phi_deriv[4];
+    double psi_deriv[4];
+    double alpha0_deriv = df_dalpha0(delta_psi, chi);
+
+    df_dtheta(thetas, phis, delta_psi, 0, theta_deriv);  // chi=0 for stability
+    df_dphi(thetas, phis, delta_psi, 0, phi_deriv);  // chi=0 for stability
+    df_dpsi(delta_psi, chi, psi_deriv);
+
+    double complex compl_velocity_conj[4];
+
+    for (int k=0; k<4; k++) {
+      // dtheta/dt
+      dy_dt[k] = -theta_deriv[k] + nu * cos(psis[k]);
+      // dphi/dt
+      dy_dt[4 + k] = -phi_deriv[k] / (sin(thetas[k]) * sin(thetas[k])) + nu * sin(psis[k]) / sin(thetas[k]);
+      // dpsi/dt
+      dy_dt[8 + k] = -1 / xi_rot * psi_deriv[k] - cos(thetas[k]) * dy_dt[4 + k];
+      // conjugated complex velocity
+      compl_velocity_conj[k] = dy_dt[k] - I * sin(thetas[k]) * dy_dt[4 + k];
+    }
+
+    // dalpha0/dt
+    double inner_summand;
+    double outer_sum=0;
+    for (int k=0; k<4; k++) {
+      inner_summand = 0;
+      for (int j=0; j<4; j++) {
+        if (j != k) {
+          inner_summand += cimag( stereogr_proj(thetas[k], phis[k]) *
+                          compl_velocity_conj[k] /
+                          (stereogr_proj(thetas[k], phis[k]) - stereogr_proj(thetas[j], phis[j])) );
+        }
+      }
+      outer_sum += 1 / 8. * ( (1 + cos(thetas[k])) * dy_dt[4 + k] + 2 / sin(thetas[k]) * inner_summand );
+    }
+    dy_dt[12] = -1 / xi_alpha * alpha0_deriv + outer_sum;
+}
+
+
+void eq_system_reduced(double* variables, double nu, double tau, double* dy_dt) {
+  // System of reduced differential equations. y=[Theta, Phi, Psi, Delta_Psi, phi1] (len(y)=5)
+
+  double theta = variables[0];
+  double phi = variables[1];
+  double avg_psi = variables[2];
+  double delta_psi = variables[3]; // - M_PI * (2 * floor((variables[3] + M_PI) / (2 * M_PI)));
+  double phi1 = variables[4];
+
+  double sin_theta = sin(theta);
+  double cos_theta = cos(theta);
+  double tan_theta = tan(theta);
+  double sin_phi = sin(phi);
+  double cos_phi = cos(phi);
+
+  // dtheta/dt
+  dy_dt[0] = nu * cos(avg_psi + delta_psi) + 2 * (4 - 8 * pow(sin_theta, 2) + 3 * pow(sin_theta, 4) * pow(sin_phi, 2)) /
+            (tan_theta * ((1 + pow(cos_theta, 2)) * (1 + pow(cos_theta, 2)) - pow(sin_theta, 4) * pow(cos_phi, 2)));
+  // dphi/dt
+  dy_dt[1] = -2 * nu * sin(avg_psi) * cos(delta_psi) / sin_theta + 8 * pow(sin_theta, 2) * sin_phi * cos_phi /
+            ((1 + pow(cos_theta, 2)) * (1 + pow(cos_theta, 2)) - pow(sin_theta, 4) * pow(cos_phi, 2));
+  // d(avg_psi)/dt
+  dy_dt[2] = cos_theta / 2 * dy_dt[1];
+  // d(delta_psi)/dt
+  dy_dt[3] = -nu * cos(avg_psi) * sin(delta_psi) / tan_theta - delta_psi / tau;
+  // dphi1/dt
+  dy_dt[4] = nu / sin_theta * cos(avg_psi) * sin(delta_psi) - dy_dt[1] / 2;
+
+}
+
+
+void propagate_trajectories_block(double* variables, double nu, double chi, double xi_rot, double xi_alpha,
+  double dt, int num_steps) {
+
+  double k1[13];
+  double k2[13];
+  double intermed[13];
+  double t = 0;
+
+  for (int i=0; i<num_steps; i++) {
+
+    eq_system(variables, nu, chi, xi_rot, xi_alpha, k1);
+
+    for (int j=0; j<13; j++) {
+      intermed[j] = variables[j] + 2. / 3. * dt * k1[j];
+    }
+
+    eq_system(intermed, nu, chi, xi_rot, xi_alpha, k2);
+
+    for (int j=0; j<13; j++) {
+      variables[j] = variables[j] + dt * (1.0 / 4.0) * (k1[j] + 3 * k2[j]);
+    }
+
+  }
+
+}
+
+
+void propagate_trajectories(double* variables0, double nu, double chi, double xi_rot, double xi_alpha,
+                            double t_span, double dt, int num_save,
+                            double* results, double* timepoints, int* nfev) {
+
+  // Should be len(results) / 13 = len(timepoints) = num_save
+
+  double t = 0;
+  int save_interval = round(t_span / (dt * (num_save - 1)));
+  dt = t_span / (save_interval * (num_save - 1));
+
+  nfev[0] = 2 * save_interval * (num_save - 1);
+
+  for (int i=0; i<num_save; i++) {
+
+    timepoints[i] = t;
+    for (int j=0; j<13; j++) {
+      results[13 * i + j] = variables0[j];
+      // printf("%f ", results[13 * i + j]);
+    }
+    // printf("\n");
+
+    if (i < num_save - 1) {
+      propagate_trajectories_block(variables0, nu, chi, xi_rot, xi_alpha, dt, save_interval);
+      t += save_interval * dt;
+    }
+
+  }
+
+}
+
+
+void propagate_reduced_trajectories_block(double* variables, double nu, double tau, double dt, int num_steps) {
+
+  double k1[5];
+  double k2[5];
+  double intermed[5];
+  double t = 0;
+
+  for (int i=0; i<num_steps; i++) {
+
+    eq_system_reduced(variables, nu, tau, k1);
+
+    for (int j=0; j<5; j++) {
+      intermed[j] = variables[j] + 2. / 3. * dt * k1[j];
+    }
+
+    eq_system_reduced(intermed, nu, tau, k2);
+
+    for (int j=0; j<5; j++) {
+      variables[j] = variables[j] + dt * (1.0 / 4.0) * (k1[j] + 3 * k2[j]);
+    }
+
+  }
+
+}
+
+
+void propagate_reduced_trajectories(double* variables0, double nu, double tau,
+                            double t_span, double dt, int num_save,
+                            double* results, double* timepoints, int* nfev) {
+
+  // Should be len(results) / 5 = len(timepoints) = num_save
+
+  double t = 0;
+  int save_interval = round(t_span / (dt * (num_save - 1)));
+  dt = t_span / (save_interval * (num_save - 1));
+
+  nfev[0] = 2 * save_interval * (num_save - 1);
+
+  for (int i=0; i<num_save; i++) {
+
+    timepoints[i] = t;
+    for (int j=0; j<5; j++) {
+      results[5 * i + j] = variables0[j];
+    }
+
+    if (i < num_save - 1) {
+      propagate_reduced_trajectories_block(variables0, nu, tau, dt, save_interval);
+      t += save_interval * dt;
+    }
+
+  }
+
+}
+
+
+void eq_system_drag(double* variables, double* grad_theta, double* grad_phi,
+                    double nu, double chi, double mu, double xi_rot, double xi_alpha, double* dy_dt) {
+    // System of differential equations. y=[theta_k, phi_k, psi_k, alpha0] (len(y)=13)
+    // grad_theta contains the etheta component of velocity gradient in radial direction
+    // grad_phi contains the ephi component of velocity gradient in radial direction
+
+
+    double thetas[4];
+    double phis[4];
+    double psis[4];
+    double delta_psi[4];
+    double alpha0 = variables[12];
+    double average_psi[4];
+
+    psi_avg(&variables[0], &variables[4], alpha0, average_psi);
+
+    for (int i=0; i<4; i++) {
+      thetas[i] = variables[i];
+      phis[i] = variables[4 + i];
+      psis[i] = variables[8 + i];
+      delta_psi[i] = psis[i] - average_psi[i];
+      delta_psi[i] -= M_PI * (2 * floor((delta_psi[i]+M_PI) / (2 * M_PI)));  // map into the interval [-pi , pi]
+    }
+
+    double theta_deriv[4];
+    double phi_deriv[4];
+    double psi_deriv[4];
+    double alpha0_deriv = df_dalpha0(delta_psi, chi);
+
+    df_dtheta(thetas, phis, delta_psi, 0, theta_deriv);  // chi=0 for stability
+    df_dphi(thetas, phis, delta_psi, 0, phi_deriv);  // chi=0 for stability
+    df_dpsi(delta_psi, chi, psi_deriv);
+
+    double complex compl_velocity_conj[4];
+
+    for (int k=0; k<4; k++) {
+      // dtheta/dt
+      dy_dt[k] = -theta_deriv[k] + nu * cos(psis[k]) - mu * grad_theta[k];
+      // dphi/dt
+      dy_dt[4 + k] = -phi_deriv[k] / (sin(thetas[k]) * sin(thetas[k])) +
+                      (nu * sin(psis[k]) - mu * grad_phi[k]) / sin(thetas[k]);
+      // dpsi/dt
+      dy_dt[8 + k] = -1 / xi_rot * psi_deriv[k] - cos(thetas[k]) * dy_dt[4 + k];
+      // conjugated complex velocity
+      compl_velocity_conj[k] = dy_dt[k] - I * sin(thetas[k]) * dy_dt[4 + k];
+    }
+
+    // dalpha0/dt
+    double inner_summand;
+    double outer_sum=0;
+    for (int k=0; k<4; k++) {
+      inner_summand = 0;
+      for (int j=0; j<4; j++) {
+        if (j != k) {
+          inner_summand += cimag( stereogr_proj(thetas[k], phis[k]) *
+                          compl_velocity_conj[k] /
+                          (stereogr_proj(thetas[k], phis[k]) - stereogr_proj(thetas[j], phis[j])) );
+        }
+      }
+      outer_sum += 1 / 8. * ( (1 + cos(thetas[k])) * dy_dt[4 + k] + 2 / sin(thetas[k]) * inner_summand );
+    }
+    dy_dt[12] = -1 / xi_alpha * alpha0_deriv + outer_sum;
+}
+
+
+void propagate_trajectories_step_drag(double* variables, double* grad_theta, double* grad_phi,
+                                      double nu, double chi, double mu, double xi_rot, double xi_alpha,
+                                      double dt) {
+  // One step of propagation with viscous drag
+
+  double k1[13];
+  double k2[13];
+  double intermed[13];
+  double t = 0;
+
+  eq_system_drag(variables, grad_theta, grad_phi, nu, chi, mu, xi_rot, xi_alpha, k1);
+
+  for (int j=0; j<13; j++) {
+    intermed[j] = variables[j] + 2. / 3. * dt * k1[j];
+  }
+
+  eq_system_drag(intermed, grad_theta, grad_phi, nu, chi, mu, xi_rot, xi_alpha, k2);
+
+  for (int j=0; j<13; j++) {
+    variables[j] = variables[j] + dt * (1.0 / 4.0) * (k1[j] + 3 * k2[j]);
+  }
+
+}
+
+
+void gaussian_belt(double* variables, double* grad_theta, double* grad_phi, double* velocity) {
+  // test function for generating the gradient field as a belt around the equator
+  double sigma = 0.3;
+  double magnitude;
+  double in_exponent;
+
+  for (int i=0; i<4; i++) {
+    in_exponent = (variables[i] - M_PI / 2) / sigma;
+    magnitude = velocity[i] / (sigma * sqrt(2 * M_PI)) * exp(-0.5 * in_exponent * in_exponent);
+    grad_theta[i] = magnitude * cos(variables[8 + i]);
+    grad_phi[i] = magnitude * sin(variables[8 + i]);
+
+  }
+
+}
+
+
+void propagate_trajectories_drag(double* variables0, double nu, double chi, double mu,
+                                double xi_rot, double xi_alpha,
+                                double t_span, double dt, int num_save, double* results,
+                                double* timepoints, int* nfev) {
+
+  // Should be len(results) / 13 = len(timepoints) = num_save
+
+  double t = 0;
+  int save_interval = round(t_span / (dt * (num_save - 1)));
+  dt = t_span / (save_interval * (num_save - 1));
+
+  double grad_theta[4];
+  double grad_phi[4];
+
+  nfev[0] = 2 * save_interval * (num_save - 1);
+
+  for (int i=0; i<num_save; i++) {
+
+    timepoints[i] = t;
+    for (int j=0; j<13; j++) {
+      results[13 * i + j] = variables0[j];
+      // printf("%f ", results[13 * i + j]);
+    }
+    // printf("\n");
+
+    double velocity[4];
+    for (int i=0; i<4; i++) {velocity[i] = nu;}
+
+    if (i < num_save - 1) {
+      for (int k=0; k<save_interval; k++) {
+        gaussian_belt(variables0, grad_theta, grad_phi, velocity);
+        propagate_trajectories_step_drag(variables0, grad_theta, grad_phi, nu, chi, mu, xi_rot, xi_alpha, dt);
+      }
+      t += save_interval * dt;
+    }
+
+  }
+
+}
+
+
+double psi_around_defect(double theta, double phi, double psi,
+  double theta_defect, double phi_defect, double psi_defect, double radial_distance) {
+  /*
+  Translation of defect velocity direction to surrounding points. Does interpolation between
+  defect direction psi_theta and previously calculated direction psi
+  */
+
+  double coord_defect[3];
+  double coord[3];
+  double axis[3];
+  double axis_normalization;
+  double cos_angle;
+  double angle;
+
+  sph_to_cart(theta_defect, phi_defect, coord_defect);
+  sph_to_cart(theta, phi, coord);
+  cross_product(coord_defect, coord, axis);
+  axis_normalization = sqrt(axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]);
+  axis[0] = axis[0] / axis_normalization;
+  axis[1] = axis[1] / axis_normalization;
+  axis[2] = axis[2] / axis_normalization;
+
+  cos_angle = cos_beta_ij(theta_defect, phi_defect, theta, phi);
+  angle = acos(cos_angle);
+
+  double defect_dir[3];
+  double vec_phi_defect[3];
+  double vec_theta_defect[3];
+  double cos_psi_defect = cos(psi_defect);
+  double sin_psi_defect = sin(psi_defect);
+
+  ephi(phi_defect, vec_phi_defect);
+  etheta(phi_defect, theta_defect, vec_theta_defect);
+
+  defect_dir[0] = cos_psi_defect * vec_theta_defect[0] + sin_psi_defect * vec_phi_defect[0];
+  defect_dir[1] = cos_psi_defect * vec_theta_defect[1] + sin_psi_defect * vec_phi_defect[1];
+  defect_dir[2] = cos_psi_defect * vec_theta_defect[2] + sin_psi_defect * vec_phi_defect[2];
+
+  // propagated (translated) defect direction
+  double propag_dir[3];
+  rodrigues_rotation(defect_dir, angle, axis, propag_dir);
+
+  // local coordiantes in he point of interest
+  double vec_phi[3];
+  double vec_theta[3];
+  double psi_propag;
+
+  ephi(phi, vec_phi);
+  etheta(phi, theta, vec_theta);
+
+  // linear interpolation if radial distance is between 0.5 and 1 in units of rho
+  // could be done with rodrigues but we used new_v = (a*v0 + (1-a)*v1) * renormalization
+  double renormalization;  // for the interpolated vector to have the same length
+  double interp_factor;
+  double cos_psi = cos(psi);
+  double sin_psi = sin(psi);
+  double previous_dir[3]; // v0 (if as v1, we consider the propagated defect direction)
+  double v0_v1;
+
+  if (radial_distance > 0.5 && radial_distance < 1) {
+    interp_factor = 2 * (1 - radial_distance);
+    previous_dir[0] = cos_psi * vec_theta[0] + sin_psi * vec_phi[0];
+    previous_dir[1] = cos_psi * vec_theta[1] + sin_psi * vec_phi[1];
+    previous_dir[2] = cos_psi * vec_theta[2] + sin_psi * vec_phi[2];
+
+    v0_v1 = previous_dir[0] * propag_dir[0] + previous_dir[1] * propag_dir[1] + previous_dir[2] * propag_dir[2];
+    renormalization = sqrt(1 / (interp_factor*interp_factor + (1 - interp_factor) * (1 - interp_factor) +
+                      2 * interp_factor * (1 - interp_factor) * v0_v1));
+
+    propag_dir[0] = (interp_factor * propag_dir[0] + (1 - interp_factor) * previous_dir[0]) * renormalization;
+    propag_dir[1] = (interp_factor * propag_dir[1] + (1 - interp_factor) * previous_dir[1]) * renormalization;
+    propag_dir[2] = (interp_factor * propag_dir[2] + (1 - interp_factor) * previous_dir[2]) * renormalization;
+
+  }
+
+  // projection to local coordinate axes
+  double sin_psi_propag = propag_dir[0] * vec_phi[0] + propag_dir[1] * vec_phi[1] + propag_dir[2] * vec_phi[2];
+  double cos_psi_propag = propag_dir[0] * vec_theta[0] + propag_dir[1] * vec_theta[1] + propag_dir[2] * vec_theta[2];
+
+  psi_propag = atan2(sin_psi_propag, cos_psi_propag);
+
+  return psi_propag;
+
+}
+
+
+void velocity_field(double* defect_config, int n, double* lattice_thetas, double* lattice_phis,
+  double nu, double rho, double* velocities) {
+
+  // printf("Slow C function.\n");
+
+  double psi_bar[4];
+  psi_avg(&defect_config[0], &defect_config[4], defect_config[12], psi_bar);
+
+  double complex z_defect[4];
+  double delta_psi[4];
+  for (int i=0; i< 4; i++) {
+    z_defect[i] = stereogr_proj(defect_config[i], defect_config[4+i]);
+    delta_psi[i] = defect_config[8 + i] - psi_bar[i];
+    delta_psi[i] = delta_psi[i] - 2 * M_PI * floor((delta_psi[i] + M_PI) / (2 * M_PI));
+  }
+
+  double alpha;
+  double dalpha_dphi, dalpha_dtheta;
+  double u_theta, u_phi;
+
+  double complex sum_in_alpha_eq;
+  double sum_in_alpha_tilde;
+  double cos_beta_defect;
+  double beta_defect;
+  double dbeta_dphi, dbeta_dtheta;
+  double alpha_eq, alpha_tilde;
+  double dalpha_eq_dphi, dalpha_eq_dtheta;
+  double dalpha_tilde_dphi, dalpha_tilde_dtheta;
+  double complex sum_in_dalpha_eq;
+  int cutoff;
+  int direction_cutoff;
+  double new_psi;  // for propagating flow direction near the defects
+  double defect_dist;  // distance from the closest defect
+  double velocity_norm;
+
+  double ln1_rho = log(1/rho);
+  double vec_phi[3];
+  double vec_theta[3];
+
+  double complex z_point;
+
+  for (int i=0; i<n; i++) {
+    sum_in_alpha_eq = 0;
+    sum_in_alpha_tilde = 0;
+    sum_in_dalpha_eq = 0;
+    cutoff = 0;
+    direction_cutoff = 0;
+    defect_dist = 0;
+
+    dalpha_tilde_dphi = 0;
+    dalpha_tilde_dtheta = 0;
+
+    // numerical correction for the lattice point exactly at the north pole (needed for the stereographic projection).
+    if (lattice_thetas[i] < 1.0E-15) {
+      lattice_thetas[i] = 1.0E-15;
+    }
+
+    z_point = stereogr_proj(lattice_thetas[i], lattice_phis[i]);
+
+    for (int j=0; j<4; j++) {
+
+      cos_beta_defect = cos_beta_ij(lattice_thetas[i], lattice_phis[i], defect_config[j], defect_config[4 + j]);
+      beta_defect = acos(cos_beta_defect);
+      if (beta_defect < 1.25 * rho) {
+        cutoff = 1;  // possible velocity cutoff near the defects (fixed defect speed, diverges otherwise)
+      }
+      if (beta_defect < rho) {
+        defect_dist = beta_defect / rho;  // distance from the defect in units of rho
+        direction_cutoff = j + 1;  // possible direction cutoff near the defects
+      }
+
+      dbeta_dphi = -dcos_beta_dphi(lattice_thetas[i], lattice_phis[i], defect_config[j], defect_config[4 + j])
+                    / sqrt(1 - cos_beta_defect * cos_beta_defect);
+      dbeta_dtheta = -dcos_beta_dtheta(lattice_thetas[i], lattice_phis[i], defect_config[j], defect_config[4 + j])
+                      / sqrt(1 - cos_beta_defect * cos_beta_defect);
+
+      sum_in_alpha_eq += clog(z_point - z_defect[j]);
+      sum_in_alpha_tilde += delta_psi[j] * (log(1 / sin(beta_defect / 2)) - 1);
+      sum_in_dalpha_eq += z_point / (z_point - z_defect[j]);
+
+      dalpha_tilde_dphi += -delta_psi[j] / (4 * ln1_rho * tan(beta_defect / 2)) * dbeta_dphi;
+      dalpha_tilde_dtheta += -delta_psi[j] / (4 * ln1_rho * tan(beta_defect / 2)) * dbeta_dtheta;
+
+    }
+
+    alpha_eq = lattice_phis[i] - defect_config[12] - 0.5 * cimag(sum_in_alpha_eq);
+    alpha_tilde = 1 / (2 * ln1_rho) * sum_in_alpha_tilde;
+    alpha = alpha_eq + alpha_tilde;
+
+    dalpha_eq_dphi = 1 - 0.5 * creal(sum_in_dalpha_eq);
+    dalpha_eq_dtheta = 1 / (2 * sin(lattice_thetas[i])) * cimag(sum_in_dalpha_eq);
+
+    dalpha_dphi = dalpha_eq_dphi + dalpha_tilde_dphi;
+    dalpha_dtheta = dalpha_eq_dtheta + dalpha_tilde_dtheta;
+
+    // averaged delta psi for higher accuracy of renormalization
+    double avg_delta_psi = 0.25 * (fabs(delta_psi[0]) + fabs(delta_psi[1]) + fabs(delta_psi[2]) + fabs(delta_psi[3]));
+
+    double normalization = 2 * rho * nu / sqrt(1 + (avg_delta_psi/ln1_rho) * (avg_delta_psi/ln1_rho));
+
+    u_theta = normalization * (sin(2 * alpha) * dalpha_dtheta
+                      - cos(2 * alpha) / sin(lattice_thetas[i]) * (cos(lattice_thetas[i]) + dalpha_dphi));
+    u_phi = normalization * (-cos(2 * alpha) * dalpha_dtheta
+                      - sin(2 * alpha) / sin(lattice_thetas[i]) * (cos(lattice_thetas[i]) + dalpha_dphi));
+
+    velocity_norm = sqrt(u_theta * u_theta + u_phi * u_phi);
+
+    // velocity cutoff near defects
+    if (cutoff == 1 && velocity_norm > nu) {
+      u_theta = u_theta * nu / velocity_norm;  // Renormalizing velocity near defects to exactly the defect speed
+      u_phi = u_phi * nu / velocity_norm;
+    }
+
+    // fixing the velocity field direction close to the defects
+    double psi_lattice;
+    if (direction_cutoff != 0) {
+      psi_lattice = atan2(u_phi, u_theta);
+
+      new_psi = psi_around_defect(lattice_thetas[i], lattice_phis[i], psi_lattice, defect_config[direction_cutoff - 1],
+        defect_config[4 + direction_cutoff - 1], defect_config[8 + direction_cutoff - 1], defect_dist);
+
+      u_theta = cos(new_psi) * nu;
+      u_phi = sin(new_psi) * nu;
+
+    }
+
+    ephi(lattice_phis[i], vec_phi);
+    etheta(lattice_phis[i], lattice_thetas[i], vec_theta);
+
+    velocities[3*i] = u_theta * vec_theta[0] + u_phi * vec_phi[0];
+    velocities[3*i+1] = u_theta * vec_theta[1] + u_phi * vec_phi[1];
+    velocities[3*i+2] = u_theta * vec_theta[2] + u_phi * vec_phi[2];
+
+  }
+
+}
+
+
+double dot_product(double* vec1,double* vec2) {
+  return vec1[0] * vec2[0] + vec1[1] * vec2[1] +vec1[2] * vec2[2];
+}
+
+
+void velocity_field_fast(double* defect_config, int n, double* lattice_thetas, double* lattice_phis,
+  double* lattice_basis, double nu, double rho, double* velocities, const int num_cores) {
+
+  // printf("Fast C function.\n");
+
+  double psi_bar[4];
+  psi_avg(&defect_config[0], &defect_config[4], defect_config[12], psi_bar);
+
+  double complex z_defect[4];
+  double delta_psi[4];
+  double defect_cart[12];
+
+  for (int i=0; i< 4; i++) {
+    z_defect[i] = stereogr_proj(defect_config[i], defect_config[4+i]);
+    delta_psi[i] = defect_config[8 + i] - psi_bar[i];
+    delta_psi[i] = delta_psi[i] - 2 * M_PI * floor((delta_psi[i] + M_PI) / (2 * M_PI));
+    sph_to_cart(defect_config[i], defect_config[4+i], &defect_cart[3*i]);
+  }
+
+  double ln1_rho = log(1/rho);
+  double vel_cutoff = cos(1.25 * rho);
+  double dir_cutoff = cos(rho);
+
+  // double start = omp_get_wtime();
+
+  #pragma omp parallel for num_threads(num_cores) schedule(dynamic, 2500)
+  for (int i=0; i<n; i++) {
+    double complex sum_in_alpha_eq = 0;
+    double complex sum_in_dalpha_eq = 0;
+    double sum_in_alpha_tilde = 0;
+    int cutoff = 0;
+    int direction_cutoff = 0;
+    double defect_dist = 0;  // distance from the closest defect
+
+    double dalpha_tilde_dphi = 0;
+    double dalpha_tilde_dtheta = 0;
+
+    double sin_lat_th = sin(lattice_thetas[i]);
+    double cos_lat_th = cos(lattice_thetas[i]);
+
+    double cos_beta_defect;
+    double sin_beta_defect;
+    double tan_beta_half;
+    double sin_beta_half;
+
+    double dbeta_dphi, dbeta_dtheta;
+
+    // numerical correction for the lattice point exactly at the north pole (needed for the stereographic projection).
+    if (lattice_thetas[i] < 1.0E-15) {
+      lattice_thetas[i] = 1.0E-15;
+    }
+
+    double complex z_point = stereogr_proj(lattice_thetas[i], lattice_phis[i]);
+
+    for (int j=0; j<4; j++) {
+
+      cos_beta_defect = dot_product(&lattice_basis[9*i], &defect_cart[3*j]);
+      if (cos_beta_defect > vel_cutoff) {
+        cutoff = 1;  // possible velocity cutoff near the defects (fixed defect speed, diverges otherwise)
+      }
+      if (cos_beta_defect > dir_cutoff) {
+        defect_dist = acos(cos_beta_defect) / rho;  // distance from the defect in units of rho
+        direction_cutoff = j + 1;  // possible direction cutoff near the defects
+      }
+
+      sin_beta_defect = sqrt(1 - cos_beta_defect * cos_beta_defect);
+
+      // derivatives over theta and phi are equal to scalar products with other basis vectors
+      dbeta_dtheta = -dot_product(&lattice_basis[9*i + 3], &defect_cart[3*j]) / sin_beta_defect;
+      dbeta_dphi = -dot_product(&lattice_basis[9*i + 6], &defect_cart[3*j]) * sin_lat_th / sin_beta_defect;
+
+      sin_beta_half = sqrt((1 - cos_beta_defect) / 2);
+      tan_beta_half = (1 - cos_beta_defect) / sin_beta_defect;
+
+      sum_in_alpha_eq += clog(z_point - z_defect[j]);
+      sum_in_alpha_tilde += delta_psi[j] * (log(1 / sin_beta_half) - 1);
+      sum_in_dalpha_eq += z_point / (z_point - z_defect[j]);
+
+      dalpha_tilde_dphi += -delta_psi[j] / (4 * ln1_rho * tan_beta_half) * dbeta_dphi;
+      dalpha_tilde_dtheta += -delta_psi[j] / (4 * ln1_rho * tan_beta_half) * dbeta_dtheta;
+
+    }
+
+    double alpha_eq = lattice_phis[i] - defect_config[12] - 0.5 * cimag(sum_in_alpha_eq);
+    double alpha_tilde = 1 / (2 * ln1_rho) * sum_in_alpha_tilde;
+    double alpha = alpha_eq + alpha_tilde;
+
+    double dalpha_eq_dphi = 1 - 0.5 * creal(sum_in_dalpha_eq);
+    double dalpha_eq_dtheta = 1 / (2 * sin_lat_th) * cimag(sum_in_dalpha_eq);
+
+    double dalpha_dphi = dalpha_eq_dphi + dalpha_tilde_dphi;
+    double dalpha_dtheta = dalpha_eq_dtheta + dalpha_tilde_dtheta;
+
+    // averaged delta psi for higher accuracy of renormalization
+    double avg_delta_psi = 0.25 * (fabs(delta_psi[0]) + fabs(delta_psi[1]) + fabs(delta_psi[2]) + fabs(delta_psi[3]));
+
+    double normalization = 2 * rho * nu / sqrt(1 + (avg_delta_psi/ln1_rho) * (avg_delta_psi/ln1_rho));
+
+    double u_theta = normalization * (sin(2 * alpha) * dalpha_dtheta
+                      - cos(2 * alpha) / sin_lat_th * (cos_lat_th + dalpha_dphi));
+    double u_phi = normalization * (-cos(2 * alpha) * dalpha_dtheta
+                      - sin(2 * alpha) / sin_lat_th * (cos_lat_th + dalpha_dphi));
+
+    double velocity_norm = sqrt(u_theta * u_theta + u_phi * u_phi);
+
+    // velocity cutoff near defects
+    if (cutoff == 1 && velocity_norm > nu) {
+      u_theta = u_theta * nu / velocity_norm;  // Renormalizing velocity near defects to exactly the defect speed
+      u_phi = u_phi * nu / velocity_norm;
+    }
+
+    // fixing the velocity field direction close to the defects
+    double psi_lattice;
+    double new_psi;  // for propagating flow direction near the defects
+    if (direction_cutoff != 0) {
+      psi_lattice = atan2(u_phi, u_theta);
+
+      new_psi = psi_around_defect(lattice_thetas[i], lattice_phis[i], psi_lattice, defect_config[direction_cutoff - 1],
+        defect_config[4 + direction_cutoff - 1], defect_config[8 + direction_cutoff - 1], defect_dist);
+
+      u_theta = cos(new_psi) * nu;
+      u_phi = sin(new_psi) * nu;
+
+    }
+
+    velocities[3*i] = u_theta * lattice_basis[9*i + 3] + u_phi * lattice_basis[9*i + 6];
+    velocities[3*i+1] = u_theta * lattice_basis[9*i + 4] + u_phi * lattice_basis[9*i + 7];
+    velocities[3*i+2] = u_theta * lattice_basis[9*i + 5] + u_phi * lattice_basis[9*i + 8];
+
+  }
+
+  // double end = omp_get_wtime();
+  // printf("Loop time: %f\n",end - start);
+
+
+}
+
+
+void velocity_field_drag(double* defect_config, int n, double* lattice_thetas, double* lattice_phis,
+  double* lattice_basis, double* lattice_grad_theta, double* lattice_grad_phi,
+  double nu, double mu, double rho, double* velocities, const int num_cores) {
+
+  // printf("Fast C function.\n");
+
+  double psi_bar[4];
+  psi_avg(&defect_config[0], &defect_config[4], defect_config[12], psi_bar);
+
+  double complex z_defect[4];
+  double delta_psi[4];
+  double defect_cart[12];
+
+  for (int i=0; i< 4; i++) {
+    z_defect[i] = stereogr_proj(defect_config[i], defect_config[4+i]);
+    delta_psi[i] = defect_config[8 + i] - psi_bar[i];
+    delta_psi[i] = delta_psi[i] - 2 * M_PI * floor((delta_psi[i] + M_PI) / (2 * M_PI));
+    sph_to_cart(defect_config[i], defect_config[4+i], &defect_cart[3*i]);
+  }
+
+  double ln1_rho = log(1/rho);
+  double vel_cutoff = cos(1.25 * rho);
+  double dir_cutoff = cos(rho);
+
+  // double start = omp_get_wtime();
+
+  #pragma omp parallel for num_threads(num_cores) schedule(dynamic, 2500)
+  for (int i=0; i<n; i++) {
+    double complex sum_in_alpha_eq = 0;
+    double complex sum_in_dalpha_eq = 0;
+    double sum_in_alpha_tilde = 0;
+    int cutoff = 0;
+    int direction_cutoff = 0;
+    double defect_dist = 0;  // distance from the closest defect
+
+    double dalpha_tilde_dphi = 0;
+    double dalpha_tilde_dtheta = 0;
+
+    double sin_lat_th = sin(lattice_thetas[i]);
+    double cos_lat_th = cos(lattice_thetas[i]);
+
+    double cos_beta_defect;
+    double sin_beta_defect;
+    double tan_beta_half;
+    double sin_beta_half;
+
+    double dbeta_dphi, dbeta_dtheta;
+
+    // numerical correction for the lattice point exactly at the north pole (needed for the stereographic projection).
+    if (lattice_thetas[i] < 1.0E-15) {
+      lattice_thetas[i] = 1.0E-15;
+    }
+
+    double complex z_point = stereogr_proj(lattice_thetas[i], lattice_phis[i]);
+
+    for (int j=0; j<4; j++) {
+
+      cos_beta_defect = dot_product(&lattice_basis[9*i], &defect_cart[3*j]);
+      if (cos_beta_defect > vel_cutoff) {
+        cutoff = 1;  // possible velocity cutoff near the defects (fixed defect speed, diverges otherwise)
+      }
+      if (cos_beta_defect > dir_cutoff) {
+        defect_dist = acos(cos_beta_defect) / rho;  // distance from the defect in units of rho
+        direction_cutoff = j + 1;  // possible direction cutoff near the defects
+      }
+
+      sin_beta_defect = sqrt(1 - cos_beta_defect * cos_beta_defect);
+
+      // derivatives over theta and phi are equal to scalar products with other basis vectors
+      dbeta_dtheta = -dot_product(&lattice_basis[9*i + 3], &defect_cart[3*j]) / sin_beta_defect;
+      dbeta_dphi = -dot_product(&lattice_basis[9*i + 6], &defect_cart[3*j]) * sin_lat_th / sin_beta_defect;
+
+      sin_beta_half = sqrt((1 - cos_beta_defect) / 2);
+      tan_beta_half = (1 - cos_beta_defect) / sin_beta_defect;
+
+      sum_in_alpha_eq += clog(z_point - z_defect[j]);
+      sum_in_alpha_tilde += delta_psi[j] * (log(1 / sin_beta_half) - 1);
+      sum_in_dalpha_eq += z_point / (z_point - z_defect[j]);
+
+      dalpha_tilde_dphi += -delta_psi[j] / (4 * ln1_rho * tan_beta_half) * dbeta_dphi;
+      dalpha_tilde_dtheta += -delta_psi[j] / (4 * ln1_rho * tan_beta_half) * dbeta_dtheta;
+
+    }
+
+    double alpha_eq = lattice_phis[i] - defect_config[12] - 0.5 * cimag(sum_in_alpha_eq);
+    double alpha_tilde = 1 / (2 * ln1_rho) * sum_in_alpha_tilde;
+    double alpha = alpha_eq + alpha_tilde;
+
+    double dalpha_eq_dphi = 1 - 0.5 * creal(sum_in_dalpha_eq);
+    double dalpha_eq_dtheta = 1 / (2 * sin_lat_th) * cimag(sum_in_dalpha_eq);
+
+    double dalpha_dphi = dalpha_eq_dphi + dalpha_tilde_dphi;
+    double dalpha_dtheta = dalpha_eq_dtheta + dalpha_tilde_dtheta;
+
+    // averaged delta psi for higher accuracy of renormalization
+    double avg_delta_psi = 0.25 * (fabs(delta_psi[0]) + fabs(delta_psi[1]) + fabs(delta_psi[2]) + fabs(delta_psi[3]));
+
+    double normalization = 2 * rho * nu / sqrt(1 + (avg_delta_psi/ln1_rho) * (avg_delta_psi/ln1_rho));
+
+    double u_theta = normalization * (sin(2 * alpha) * dalpha_dtheta
+                      - cos(2 * alpha) / sin_lat_th * (cos_lat_th + dalpha_dphi));
+    double u_phi = normalization * (-cos(2 * alpha) * dalpha_dtheta
+                      - sin(2 * alpha) / sin_lat_th * (cos_lat_th + dalpha_dphi));
+
+    double velocity_norm = sqrt(u_theta * u_theta + u_phi * u_phi);
+
+    // velocity cutoff near defects
+    if (cutoff == 1 && velocity_norm > nu) {
+      u_theta = u_theta * nu / velocity_norm;  // Renormalizing velocity near defects to exactly the defect speed
+      u_phi = u_phi * nu / velocity_norm;
+    }
+
+    // fixing the velocity field direction close to the defects
+    double psi_lattice;
+    double new_psi;  // for propagating flow direction near the defects
+    if (direction_cutoff != 0) {
+      psi_lattice = atan2(u_phi, u_theta);
+
+      new_psi = psi_around_defect(lattice_thetas[i], lattice_phis[i], psi_lattice, defect_config[direction_cutoff - 1],
+        defect_config[4 + direction_cutoff - 1], defect_config[8 + direction_cutoff - 1], defect_dist);
+
+      u_theta = cos(new_psi) * nu;
+      u_phi = sin(new_psi) * nu;
+
+    }
+
+    u_theta = u_theta - mu * lattice_grad_theta[i];
+    u_phi = u_phi - mu * lattice_grad_phi[i];
+
+    velocities[3*i] = u_theta * lattice_basis[9*i + 3] + u_phi * lattice_basis[9*i + 6];
+    velocities[3*i+1] = u_theta * lattice_basis[9*i + 4] + u_phi * lattice_basis[9*i + 7];
+    velocities[3*i+2] = u_theta * lattice_basis[9*i + 5] + u_phi * lattice_basis[9*i + 8];
+
+  }
+
+  // double end = omp_get_wtime();
+  // printf("Loop time: %f\n",end - start);
+
+
+}

dipol/fibonacci/saturn_ring/state_A1_test/konstante.c

@@ -0,0 +1,82 @@
+//Logical markers for the datapoint characterisation
+#define LMARK1			2				//bulk
+#define LMARKKOL		4				//colloid
+#define LMARKPOV		8				//surface node
+#define LMARKRP			16				//Marker is a part of boundary conditions
+#define LMARK0			32				//There is absolutely nothing to be done in that point
+
+
+//Numerical parameters
+
+//A time step in the code corresponds to RESP^2/L/GAM*DT, where RESP=resolution, L=elastic constant, GAM is GAMMA viscosity. All can be chosen arbitrarily
+#define DT			0.025  //Changing DT does not alter material parameter values (but it affects the Reynolds number)
+
+//relaxation time of the Lb scheme -- sets up the isotropic visosity alpha4. It is important for numerical precision, in particular the slip velocity at the walls.
+#define TAUF			(2.f*DT)
+
+//Density parameter - vaguely related to the actual density (note that LB is weakly compressible). Sets up the isotropic viscosity eta. The fact that the density is way too largre results in larger Reynolds number compared to experiments. eta=DENSITYINIT*DT*(TAUF/DT-0.5)/3 in units of 1/GAMMA. Lower DENSITYINIT might lead to numerical instabilities.
+#define	DENSITYINIT		(2.75f/DT)
+
+
+
+//Material parameters
+
+//Alignment parameter
+#define XI			1.f				//Tumbling/aligning regime -- transition is at XI=0.857 (at SEQ=0.533)
+
+//Activity
+float ZETA;//			0.0f				//Activity parameter. To get the actual value (in units of L/RESP^2) multiply by RESP^2/L
+
+//Phase parameters (in the units of L/RESP^2)
+#define A			(-0.43f)
+#define B			(-5.3f)
+#define C			(4.325f)
+#define SEQ			(0.5*(-B/(3*C)+sqrt(fsqr(B/(3*C))-(8*A)/(3*C))) )
+//For the above phase parameters, nematic correlation length KSI=sqrt(L/(A+B*Seq+9/2*C*Seq^2)) equals (0.663e-8) times RESP
+
+// Surface anchoring (in the units of L/RESP)
+#define WHOM			(0.5f)			//homeotropic anchoring strength
+// #define WHOM			(0.005f)			//homeotropic anchoring strength
+
+
+// Radius of a shell
+#define R			(30.0f)			
+
+//active droplet
+#define dt_active0  (1.0E-4)  // ACTIVE DROPLET TIME/PASSIVE NEMATIC TIME  
+#define VELOCITY_SCALING  (0.01)  
+
+
+// 
+// #define nu			  5
+// #define rho			(0.1f)
+// #define chi			1.0/log(1.0/rho)
+// #define xi_rot			1.0
+// #define tau_active			xi_rot / (2 * chi)
+// #define xi_alpha			1.0
+
+
+
+
+//A possible choice for the material parameters can correspond to the values that Miha Ravnik chose long time ago
+// //parameters that define the scale of the system (not really, these are just a suggestion)
+// #define RESP			(1e-8)				//Resolution (in meters)
+// //#define KSI			(0.663e-8)
+// #define L			(4.e-11)			//Elastic constant -- (in N)
+// //Phase parameters
+// #define SEQ			(0.5*(-B/(3*C)+sqrt(fsqr(B/(3*C))-(8*A)/(3*C))) )
+// #define Ape			(-0.172e6)			// J/m3
+// #define Bpe			(-2.12e6)			// J/m3
+// #define Cpe			(1.73e6)			// J/m3
+// #define A			(Ape*RESP*RESP/L)		// (for T-T*=-10K)   A=-0.43
+// #define B			(Bpe*RESP*RESP/L)		// -5.3
+// #define C			(Cpe*RESP*RESP/L)		// 4.325
+
+
+
+
+
+//Pressure flow driving (typically not needed)
+// #define DELTAP			0
+// #define DELTAPX			(0.001f*DENSITYINIT)
+// #define DELTAPY			(0.0f*DENSITYINIT)

dipol/fibonacci/saturn_ring/state_A1_test/num_rec.c

@@ -0,0 +1,589 @@
+//Useful mathematical functions from numerical recipes library
+
+//Square function for floats and doubles
+#ifndef SQR_H
+#define SQR_H
+
+#ifdef __GNUC__
+#define SQR_ATTR __attribute__((__const__))
+#else
+#define SQR_ATTR
+#endif
+
+
+SQR_ATTR double sqr(double a)
+  {
+   return a*a;
+  }
+
+SQR_ATTR float fsqr(float a)
+  {
+   return a*a;
+  }
+
+
+#undef SQR_ATTR
+
+#endif
+
+
+#ifndef _NR_UTILS_H_
+#define _NR_UTILS_H_
+
+static float sqrarg;
+#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
+
+static double dsqrarg;
+#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
+
+static double dmaxarg1,dmaxarg2;
+#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
+        (dmaxarg1) : (dmaxarg2))
+
+static double dminarg1,dminarg2;
+#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\
+        (dminarg1) : (dminarg2))
+
+static float maxarg1,maxarg2;
+#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
+        (maxarg1) : (maxarg2))
+
+static float minarg1,minarg2;
+#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\
+        (minarg1) : (minarg2))
+
+static long lmaxarg1,lmaxarg2;
+#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\
+        (lmaxarg1) : (lmaxarg2))
+
+static long lminarg1,lminarg2;
+#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\
+        (lminarg1) : (lminarg2))
+
+static int imaxarg1,imaxarg2;
+#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\
+        (imaxarg1) : (imaxarg2))
+
+static int iminarg1,iminarg2;
+#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
+        (iminarg1) : (iminarg2))
+
+#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
+
+void nrerror(char error_text[]);
+float *vector(long nl, long nh);
+int *ivector(long nl, long nh);
+unsigned char *cvector(long nl, long nh);
+unsigned long *lvector(long nl, long nh);
+double *dvector(long nl, long nh);
+float **matrix(long nrl, long nrh, long ncl, long nch);
+double **dmatrix(long nrl, long nrh, long ncl, long nch);
+int **imatrix(long nrl, long nrh, long ncl, long nch);
+float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
+	long newrl, long newcl);
+float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch);
+float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
+void free_vector(float *v, long nl, long nh);
+void free_ivector(int *v, long nl, long nh);
+void free_cvector(unsigned char *v, long nl, long nh);
+void free_lvector(unsigned long *v, long nl, long nh);
+void free_dvector(double *v, long nl, long nh);
+void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
+void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
+void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
+void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch);
+void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch);
+void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
+	long ndl, long ndh);
+
+#endif /* _NR_UTILS_H_ */
+#include <stdio.h>
+#include <stddef.h>
+#include <stdlib.h>
+#define NR_END 1
+#define FREE_ARG char*
+
+void nrerror(char error_text[])
+/* Numerical Recipes standard error handler */
+{
+	fprintf(stderr,"Numerical Recipes run-time error...\n");
+	fprintf(stderr,"%s\n",error_text);
+	fprintf(stderr,"...now exiting to system...\n");
+	exit(1);
+}
+
+float *vector(long nl, long nh)
+/* allocate a float vector with subscript range v[nl..nh] */
+{
+	float *v;
+
+	v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
+	if (!v) nrerror("allocation failure in vector()");
+	return v-nl+NR_END;
+}
+
+int *ivector(long nl, long nh)
+/* allocate an int vector with subscript range v[nl..nh] */
+{
+	int *v;
+
+	v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
+	if (!v) nrerror("allocation failure in ivector()");
+	return v-nl+NR_END;
+}
+
+unsigned char *cvector(long nl, long nh)
+/* allocate an unsigned char vector with subscript range v[nl..nh] */
+{
+	unsigned char *v;
+
+	v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
+	if (!v) nrerror("allocation failure in cvector()");
+	return v-nl+NR_END;
+}
+
+unsigned long *lvector(long nl, long nh)
+/* allocate an unsigned long vector with subscript range v[nl..nh] */
+{
+	unsigned long *v;
+
+	v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long)));
+	if (!v) nrerror("allocation failure in lvector()");
+	return v-nl+NR_END;
+}
+
+double *dvector(long nl, long nh)
+/* allocate a double vector with subscript range v[nl..nh] */
+{
+	double *v;
+
+	v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
+	if (!v) nrerror("allocation failure in dvector()");
+	return v-nl+NR_END;
+}
+
+float **matrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	float **m;
+
+	/* allocate pointers to rows */
+	m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+double **dmatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	double **m;
+
+	/* allocate pointers to rows */
+	m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+int **imatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	int **m;
+
+	/* allocate pointers to rows */
+	m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
+	long newrl, long newcl)
+/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
+{
+	long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
+	float **m;
+
+	/* allocate array of pointers to rows */
+	m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
+	if (!m) nrerror("allocation failure in submatrix()");
+	m += NR_END;
+	m -= newrl;
+
+	/* set pointers to rows */
+	for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
+/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
+declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
+and ncol=nch-ncl+1. The routine should be called with the address
+&a[0][0] as the first argument. */
+{
+	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	float **m;
+
+	/* allocate pointers to rows */
+	m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
+	if (!m) nrerror("allocation failure in convert_matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* set pointers to rows */
+	m[nrl]=a-ncl;
+	for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
+/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
+{
+	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
+	float ***t;
+
+	/* allocate pointers to pointers to rows */
+	t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
+	if (!t) nrerror("allocation failure 1 in f3tensor()");
+	t += NR_END;
+	t -= nrl;
+
+	/* allocate pointers to rows and set pointers to them */
+	t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
+	if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
+	t[nrl] += NR_END;
+	t[nrl] -= ncl;
+
+	/* allocate rows and set pointers to them */
+	t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
+	if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
+	t[nrl][ncl] += NR_END;
+	t[nrl][ncl] -= ndl;
+
+	for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
+	for(i=nrl+1;i<=nrh;i++) {
+		t[i]=t[i-1]+ncol;
+		t[i][ncl]=t[i-1][ncl]+ncol*ndep;
+		for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
+	}
+
+	/* return pointer to array of pointers to rows */
+	return t;
+}
+
+void free_vector(float *v, long nl, long nh)
+/* free a float vector allocated with vector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_ivector(int *v, long nl, long nh)
+/* free an int vector allocated with ivector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_cvector(unsigned char *v, long nl, long nh)
+/* free an unsigned char vector allocated with cvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_lvector(unsigned long *v, long nl, long nh)
+/* free an unsigned long vector allocated with lvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_dvector(double *v, long nl, long nh)
+/* free a double vector allocated with dvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
+/* free a float matrix allocated by matrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
+/* free a double matrix allocated by dmatrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
+/* free an int matrix allocated by imatrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
+/* free a submatrix allocated by submatrix() */
+{
+	free((FREE_ARG) (b+nrl-NR_END));
+}
+
+void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
+/* free a matrix allocated by convert_matrix() */
+{
+	free((FREE_ARG) (b+nrl-NR_END));
+}
+
+void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
+	long ndl, long ndh)
+/* free a float f3tensor allocated by f3tensor() */
+{
+	free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
+	free((FREE_ARG) (t[nrl]+ncl-NR_END));
+	free((FREE_ARG) (t+nrl-NR_END));
+}
+#include <math.h>
+#define NRANSI
+
+
+float pythag(float a, float b)
+{
+	float absa,absb;
+	absa=fabs(a);
+	absb=fabs(b);
+	if (absa > absb) return absa*sqrt(1.0+SQR(absb/absa));
+	else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb)));
+}
+#undef NRANSI
+/* note #undef's at end of file */
+#define IM1 2147483563
+#define IM2 2147483399
+#define AM (1.0/IM1)
+#define IMM1 (IM1-1)
+#define IA1 40014
+#define IA2 40692
+#define IQ1 53668
+#define IQ2 52774
+#define IR1 12211
+#define IR2 3791
+#define NTAB 32
+#define NDIV (1+IMM1/NTAB)
+#define EPS 1.2e-7
+#define RNMX (1.0-EPS)
+
+float ran2(long *idum)
+{
+	int j;
+	long k;
+	static long idum2=123456789;
+	static long iy=0;
+	static long iv[NTAB];
+	float temp;
+
+	if (*idum <= 0) {
+		if (-(*idum) < 1) *idum=1;
+		else *idum = -(*idum);
+		idum2=(*idum);
+		for (j=NTAB+7;j>=0;j--) {
+			k=(*idum)/IQ1;
+			*idum=IA1*(*idum-k*IQ1)-k*IR1;
+			if (*idum < 0) *idum += IM1;
+			if (j < NTAB) iv[j] = *idum;
+		}
+		iy=iv[0];
+	}
+	k=(*idum)/IQ1;
+	*idum=IA1*(*idum-k*IQ1)-k*IR1;
+	if (*idum < 0) *idum += IM1;
+	k=idum2/IQ2;
+	idum2=IA2*(idum2-k*IQ2)-k*IR2;
+	if (idum2 < 0) idum2 += IM2;
+	j=iy/NDIV;
+	iy=iv[j]-idum2;
+	iv[j] = *idum;
+	if (iy < 1) iy += IMM1;
+	if ((temp=AM*iy) > RNMX) return RNMX;
+	else return temp;
+}
+#undef IM1
+#undef IM2
+#undef AM
+#undef IMM1
+#undef IA1
+#undef IA2
+#undef IQ1
+#undef IQ2
+#undef IR1
+#undef IR2
+#undef NTAB
+#undef NDIV
+#undef EPS
+#undef RNMX
+#include <math.h>
+#define NRANSI
+
+int tqli(float d[], float e[], int n, float **z)
+{
+	float pythag(float a, float b);
+	int m,l,iter,i,k;
+	float s,r,p,g,f,dd,c,b;
+
+	for (i=2;i<=n;i++) e[i-1]=e[i];
+	e[n]=0.0;
+	for (l=1;l<=n;l++) {
+		iter=0;
+		do {
+			for (m=l;m<=n-1;m++) {
+				dd=fabs(d[m])+fabs(d[m+1]);
+				if ((float)(fabs(e[m])+dd) == dd) break;
+			}
+			if (m != l) {
+				if (iter++ == 30){
+                    nrerror("Too many iterations in tqli");
+                    return 1;
+                }
+				g=(d[l+1]-d[l])/(2.0*e[l]);
+				r=pythag(g,1.0);
+				g=d[m]-d[l]+e[l]/(g+SIGN(r,g));
+				s=c=1.0;
+				p=0.0;
+				for (i=m-1;i>=l;i--) {
+					f=s*e[i];
+					b=c*e[i];
+					e[i+1]=(r=pythag(f,g));
+					if (r == 0.0) {
+						d[i+1] -= p;
+						e[m]=0.0;
+						break;
+					}
+					s=f/r;
+					c=g/r;
+					g=d[i+1]-p;
+					r=(d[i]-g)*s+2.0*c*b;
+					d[i+1]=g+(p=s*r);
+					g=c*r-b;
+					for (k=1;k<=n;k++) {
+						f=z[k][i+1];
+						z[k][i+1]=s*z[k][i]+c*f;
+						z[k][i]=c*z[k][i]-s*f;
+					}
+				}
+				if (r == 0.0 && i >= l) continue;
+				d[l] -= p;
+				e[l]=g;
+				e[m]=0.0;
+			}
+		} while (m != l);
+	}
+	return 0;
+}
+#undef NRANSI
+#include <math.h>
+
+void tred2(float **a, int n, float d[], float e[])
+{
+	int l,k,j,i;
+	float scale,hh,h,g,f;
+
+	for (i=n;i>=2;i--) {
+		l=i-1;
+		h=scale=0.0;
+		if (l > 1) {
+			for (k=1;k<=l;k++)
+				scale += fabs(a[i][k]);
+			if (scale == 0.0)
+				e[i]=a[i][l];
+			else {
+				for (k=1;k<=l;k++) {
+					a[i][k] /= scale;
+					h += a[i][k]*a[i][k];
+				}
+				f=a[i][l];
+				g=(f >= 0.0 ? -sqrt(h) : sqrt(h));
+				e[i]=scale*g;
+				h -= f*g;
+				a[i][l]=f-g;
+				f=0.0;
+				for (j=1;j<=l;j++) {
+					a[j][i]=a[i][j]/h;
+					g=0.0;
+					for (k=1;k<=j;k++)
+						g += a[j][k]*a[i][k];
+					for (k=j+1;k<=l;k++)
+						g += a[k][j]*a[i][k];
+					e[j]=g/h;
+					f += e[j]*a[i][j];
+				}
+				hh=f/(h+h);
+				for (j=1;j<=l;j++) {
+					f=a[i][j];
+					e[j]=g=e[j]-hh*f;
+					for (k=1;k<=j;k++)
+						a[j][k] -= (f*e[k]+g*a[i][k]);
+				}
+			}
+		} else
+			e[i]=a[i][l];
+		d[i]=h;
+	}
+	d[1]=0.0;
+	e[1]=0.0;
+	/* Contents of this loop can be omitted if eigenvectors not
+			wanted except for statement d[i]=a[i][i]; */
+	for (i=1;i<=n;i++) {
+		l=i-1;
+		if (d[i]) {
+			for (j=1;j<=l;j++) {
+				g=0.0;
+				for (k=1;k<=l;k++)
+					g += a[i][k]*a[k][j];
+				for (k=1;k<=l;k++)
+					a[k][j] -= g*a[k][i];
+			}
+		}
+		d[i]=a[i][i];
+		a[i][i]=1.0;
+		for (j=1;j<=l;j++) a[j][i]=a[i][j]=0.0;
+	}
+}