1 files changed, 167 insertions, 0 deletions
diff --git a/src/flow.c b/src/flow.c
new file mode 100644
index 0000000..d271e44
--- /dev/null
+++ b/src/flow.c
@@ -0,0 +1,167 @@
+/*
+Copyright 2015 Ian Jauslin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "flow.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include "tools.h"
+#include "math.h"
+#include "definitions.cpp"
+#include "number.h"
+#include "array.h"
+#include "coefficient.h"
+
+
+
+// compute flow numerically, no exponentials
+//    inputs: flow_equation
+//            init, niter, tol (the allowed error at each step), ls (whether to display the results in terms of ls), display_mode (what to print)
+int numerical_flow(Grouped_Polynomial flow_equation, RCC init, Labels labels, int niter, long double tol, int display_mode){
+  // running coupling contants
+  RCC rccs=init;
+  int i,j;
+
+  if(display_mode==DISPLAY_NUMERICAL){
+    // print labels
+    printf("%5s  ","n");
+    for(j=0;j<rccs.length;j++){
+      print_label(rccs.indices[j], labels);
+    }
+    printf("\n\n");
+
+    // print initial values
+    printf("%5d  ",0);
+    for(j=0;j<rccs.length;j++){
+      printf("% 14.7Le  ",rccs.values[j]);
+    }
+    printf("\n");
+  }
+
+  for(i=0;i<niter;i++){
+    // compute a single step
+    step_flow(&rccs, flow_equation, tol);
+    // convert ls to alphas
+    if(display_mode==DISPLAY_NUMERICAL){
+      // print the result
+      printf("%5d  ",i+1);
+      for(j=0;j<rccs.length;j++){
+	printf("% 14.7Le  ",rccs.values[j]);
+      }
+      printf("\n");
+    }
+  }
+
+  if(display_mode==DISPLAY_NUMERICAL){
+    // print labels
+    printf("\n");
+    printf("%5s  ","n");
+    for(j=0;j<rccs.length;j++){
+      print_label(rccs.indices[j], labels);
+    }
+    printf("\n\n");
+  }
+
+  if(display_mode==DISPLAY_FINAL){
+    RCC_print(rccs);
+  }
+
+  return(0);
+}
+
+// single step in the flow no exponentials
+//          inputs: flow_equation, tol
+//   input/outputs: rccs
+int step_flow(RCC* rccs, Grouped_Polynomial flow_equation, long double tol){
+  int i;
+  long double* new_rccs=calloc((*rccs).length,sizeof(long double));
+  Int_Array computed;
+
+  init_Int_Array(&computed, (*rccs).length);
+
+  // initialize vectors to 0
+  for(i=0;i<(*rccs).length;i++){
+    new_rccs[i]=0.;
+  }
+
+  // compute the constants first
+  for(i=0;i<flow_equation.length;i++){
+    if(flow_equation.indices[i]<0){
+      evalcoef(*rccs, flow_equation.coefs[i], new_rccs+i);
+      // if the new rcc is too small, then ignore it
+      if(fabs(new_rccs[i])<tol){
+	new_rccs[i]=0.;
+      }
+      (*rccs).values[i]=new_rccs[i];
+    }
+  }
+
+  // for each equation
+  for(i=0;i<flow_equation.length;i++){
+    if(flow_equation.indices[i]>=0){
+      evalcoef(*rccs, flow_equation.coefs[i], new_rccs+i);
+      // if the new rcc is too small, then ignore it
+      if(fabs(new_rccs[i])<tol){
+	new_rccs[i]=0.;
+      }
+    }
+  }
+
+  // copy results to rccs
+  for(i=0;i<(*rccs).length;i++){
+    (*rccs).values[i]=new_rccs[i];
+  }
+
+  // free memory
+  free_Int_Array(computed);
+  free(new_rccs);
+  return(0);
+}
+
+
+// print the label of an rcc (takes constants and derivatives into account)
+int print_label(int index, Labels labels){
+  int i;
+  int nderivs;
+  int posin_labels;
+  Char_Array label;
+
+  // constant term
+  if(index<0){
+    nderivs=-index/DOFFSET;
+    for (i=0;i<12-nderivs;i++){
+      printf(" ");
+    }
+    for(i=0;i<nderivs;i++){
+      printf("d");
+    }
+    printf("C%d  ",-index-nderivs*DOFFSET);
+  }
+  else{
+    nderivs=index/DOFFSET;
+    posin_labels=intlist_find_err(labels.indices, labels.length, index-nderivs*DOFFSET);
+    label=labels.labels[posin_labels];
+    for (i=0;i<14-label.length-nderivs;i++){
+      printf(" ");
+    }
+    for(i=0;i<nderivs;i++){
+      printf("d");
+    }
+    printf("%s  ",char_array_to_str_noinit(labels.labels+posin_labels));
+  }
+
+  return(0);
+}