diff --git a/make_dyadic_vectors.cc b/make_dyadic_vectors.cc
index aee608a0c8ece5dccea9814d78c61a2e626ae4f3..3d81874d568c69b9d68c46d5081d53d4967e07ff 100644
--- a/make_dyadic_vectors.cc
+++ b/make_dyadic_vectors.cc
@@ -9,10 +9,12 @@
 using namespace std;
 using namespace NEWIMAGE;
 
+
 int main ( int argc, char **argv ){
   if(argc<4){
-    cerr<<"usage: make_dyadic_vectors <theta_vol> <phi_vol> [mask] <output>"<<endl;
-    cerr<<"[mask] is optional"<<endl;
+    cout<<"usage: make_dyadic_vectors <theta_vol> <phi_vol> [mask] <output> [perc]"<<endl;
+    cout<<"[mask] is optional"<<endl;
+    cout<<"[perc] is the {perc}% angle of the output cone of uncertainty (output will be in degrees)"<<endl;
     exit(1);
   }
 
@@ -21,17 +23,34 @@ int main ( int argc, char **argv ){
   read_volume4D(phs,argv[2]);
   volume<float> mask;
   string oname;
-  if(argc==5){
-    oname=argv[4];
-    read_volume(mask,argv[3]);
-  }
-  else{
+  bool ocone  = false;
+  float acone = 0;
+  if(argc==4){
     mask=ths[0]*0+1;
     oname=argv[3];
   }
+  else if(argc>=5){
+    oname=argv[4];
+    read_volume(mask,argv[3]);
+    if(argc==6){
+      ocone = true;
+      acone = atof(argv[5]);
+      acone/=100.0;
+      if(acone<0 || acone>1){
+	cerr << "cone of uncertainty should be in percent and range between 0 and 100%" << endl;
+	exit(1);
+      }
+
+    }
+  }
+
+  
+
   volume4D<float> dyadic_vecs(ths.xsize(),ths.ysize(),ths.zsize(),3);
+  volume<float> disp(ths.xsize(),ths.ysize(),ths.zsize());
   dyadic_vecs=0;
   copybasicproperties(ths[0],dyadic_vecs[0]);
+  copybasicproperties(ths[0],disp);
   SymmetricMatrix dyad(3);dyad=0;
   ColumnVector dir(3);
   
@@ -41,7 +60,7 @@ int main ( int argc, char **argv ){
   for(int k=ths.minz();k<=ths.maxz();k++){
     for(int j=ths.miny();j<=ths.maxy();j++){
       for(int i=ths.minx();i<=ths.maxx();i++){
-	if(mask(i,j,k)>0){
+	if(mask(i,j,k)!=0){
 	  dyad=0;
 	  for(int s=ths.mint();s<=ths.maxt();s++){		   
 	  
@@ -52,7 +71,7 @@ int main ( int argc, char **argv ){
 	    dir(3)=cos(th);
 	    dyad << dyad+dir*dir.t();
 	  }
-	  
+	  dyad = dyad/float(ths.maxt()-ths.mint()+1);
 	  EigenValues(dyad,dyad_D,dyad_V);
 	  int maxeig;
 	  if(dyad_D(1)>dyad_D(2)){
@@ -66,17 +85,53 @@ int main ( int argc, char **argv ){
 	  dyadic_vecs(i,j,k,0)=dyad_V(1,maxeig);
 	  dyadic_vecs(i,j,k,1)=dyad_V(2,maxeig);
 	  dyadic_vecs(i,j,k,2)=dyad_V(3,maxeig);
+	  disp(i,j,k)=1-dyad_D.MaximumAbsoluteValue();
 	}
 	else{
 	  dyadic_vecs(i,j,k,0)=0;
 	  dyadic_vecs(i,j,k,1)=0;
 	  dyadic_vecs(i,j,k,2)=0;
+	  disp(i,j,k)=0;
 	}
       }
     }
   }
-  
+
   save_volume4D(dyadic_vecs,oname);
+  save_volume(disp,oname+"_dispersion");
+
+  // where we calculate the cone of uncertainty
+  if(ocone){
+    cout << "calculate cones" << endl;
+    volume<float> cones;
+    cones = ths[0]*0;
+    ColumnVector meanDyad(3);
+    for(int k=ths.minz();k<=ths.maxz();k++)
+      for(int j=ths.miny();j<=ths.maxy();j++)
+	for(int i=ths.minx();i<=ths.maxx();i++){
+	  if(mask(i,j,k)==0)continue;
+	  
+	  meanDyad<< dyadic_vecs(i,j,k,0) << dyadic_vecs(i,j,k,1) << dyadic_vecs(i,j,k,2);
+	  ColumnVector angles(ths.tsize());
+	  for(int s=0;s<ths.tsize();s++){		   
+	  
+	    float th=ths(i,j,k,s);
+	    float ph=phs(i,j,k,s);
+	    dir(1)=sin(th)*cos(ph);
+	    dir(2)=sin(th)*sin(ph);
+	    dir(3)=cos(th);
+
+	    angles(s+1) = acos( abs(dot(dir,meanDyad)) ) * 180.0/M_PI;
+	  }
+	  SortAscending(angles);
+	  cones(i,j,k) = angles( round(angles.Nrows()*acone) );
+
+	}
+  
+    save_volume(cones,oname+"_cones"+num2str(acone*100));
+  }
+
+  return 0;
 }