Permalink
Browse files

Increased default precision of pcgtol, the conjugate gradient parameter.

This is necessary for larger systems.

Changed the turbulence generation file to scale turbulence correctly when there is an outer scale.

A number of minor changes to yao.i
  • Loading branch information...
1 parent 56d5239 commit 267aae87fd5ce4f95ec593ab435fd1221638d336 Marcos committed Nov 1, 2011
Showing with 26 additions and 21 deletions.
  1. +1 −1 aoutil.i
  2. +15 −7 turbulence.i
  3. +10 −13 yao.i
View
@@ -654,7 +654,7 @@ func check_parameters(void)
if (mat.sparse_MR == long()){mat.sparse_MR = 10000;}
if (mat.sparse_MN == long()){mat.sparse_MN = 200000;}
if (mat.sparse_thresh == float()){mat.sparse_thresh = 1e-8;}
- if (mat.sparse_pcgtol == float()){mat.sparse_pcgtol = 1e-3;}
+ if (mat.sparse_pcgtol == float()){mat.sparse_pcgtol = 1e-6;}
if (mat.fit_subsamp == long()){mat.fit_subsamp = 1;}
if (mat.fit_target == long()){mat.fit_target = 1;}
if (mat.fit_minval == float()){mat.fit_minval = 1e-2;}
View
@@ -116,18 +116,26 @@ func create_phase_screens(dimx,dimy,l0=,prefix=,nalias=,no_ipart=)
psfunc(i) = sqrt((fsx+fsy)/2.);
gui_progressbar_frac,0.25+0.6*(i-off(1))/(off(2)-off(1));
}
- theo = sqrt(6.88*indgen(off(2))^1.66);
-
- // write,format="normalization factor = %f\n",sqrt(6.88*off^1.666/fs);
+
+ c = (24./5.*gamma(6/5.))^(5/6.);
+ r = float(indgen(off(2)));
+ if (l0 == 0){
+ theo = sqrt(2*c*r^(5./3.));
+ } else {
+ f0 = 1./l0;
+ c = (24./5.*gamma(6/5.))^(5/6.);
+ r = float(indgen(off(2)));
+ theo = sqrt(2*c*gamma(11./6.)/(2^(5./6)*pi^(8./3))*(f0)^(-5./3)*(gamma(5./6.)/2^(1./6.) - (2*pi*r*f0)^(5./6.)*gsl_sf_bessel_Knu(5./6., 2*pi*r*f0)));
+ }
+
write,format="normalization factor (actual/theo)= %f\n",
(nfact=avg(psfunc(off(1):off(2))/theo(off(1):off(2))));
write,psfunc(off(1):off(2))/theo(off(1):off(2));
pscreen(*) = pscreen(*)/float(nfact);
-
- print,"Sectionning and saving phase screens";
- gui_progressbar_text,"Sectionning and saving phase screens";
+ print,"Sectioning and saving phase screens";
+ gui_progressbar_text,"Sectioning and saving phase screens";
pscreen = reform(pscreen,[3,dimx,dimy,nscreen]);
if (!is_void(prefix)) {
for (i=1;i<=nscreen;i++) {
@@ -188,7 +196,7 @@ func generate_phase(dim)
func generate_von_karman_spectrum(dim,k0,nalias=,silent=)
/* DOCUMENT func generate_von_karman_spectrum(sdim,bdim,k0)
- generate correct VoKarman spectrum including aliasing.
+ generate correct von Karman spectrum including aliasing.
SEE ALSO:
*/
View
@@ -883,7 +883,7 @@ func build_cmat(condition,modalgain,subsystem=,all=,nomodalgain=,disp=)
F.Rigaut, June 17,2002
condition = Filter eigenvalues ev (and modes) for max(ev)/ev > condition
- modalgain = vector of system mode gains to pre-multiply the contolr matrix
+ modalgain = vector of system mode gains to pre-multiply the control matrix
with. rarely used as these modes are not generally natural
w.r.t. the turbulence.
all = if not set, one (1) mode is forced to be discarded (useful if regular
@@ -2425,7 +2425,7 @@ func aoinit(disp=,clean=,forcemat=,svd=,dpi=,keepdmconfig=)
} while (again != "n");
}
- if (dm(nm).fit_wfs){// if DM on WFS path only
+ if (dm(nm).fit_wfs){// if DM used for WFS
ok = indgen(numberof(tmp));
nok = [];
} else {
@@ -2749,7 +2749,7 @@ func aoinit(disp=,clean=,forcemat=,svd=,dpi=,keepdmconfig=)
for (nm=1;nm<=numberof(dm);nm++){
if (dm(nm).fitvirtualdm){
- virtualDMs = *dm(nm).fitvirtualdm;
+ virtualDMs = int(*dm(nm).fitvirtualdm);
nVirtualDMs = numberof(virtualDMs);
if (dm(nm).type == "stackarray"){
@@ -2803,7 +2803,7 @@ func aoinit(disp=,clean=,forcemat=,svd=,dpi=,keepdmconfig=)
for (nm=1;nm<=numberof(dm);nm++){
if (dm(nm).fitvirtualdm){
- virtualDMs = *dm(nm).fitvirtualdm;
+ virtualDMs = int(*dm(nm).fitvirtualdm);
nVirtualDMs = numberof(virtualDMs);
// calculate the matrix that shows how the tomographic DM affects the phase
@@ -2996,7 +2996,7 @@ func aoinit(disp=,clean=,forcemat=,svd=,dpi=,keepdmconfig=)
} else {
// tomographic DM
vidx = [];
- virtualDMs = *dm(nm).fitvirtualdm;
+ virtualDMs = int(*dm(nm).fitvirtualdm);
for (c1=1;c1<= numberof(virtualDMs);c1++){
grow, vidx, indgen(indexDm(1,virtualDMs(c1)):indexDm(2,virtualDMs(c1)));
}
@@ -3134,8 +3134,7 @@ func aoinit(disp=,clean=,forcemat=,svd=,dpi=,keepdmconfig=)
}
}
- fMatSP = rcotr(fMatSP);
- t1 = rcoatb(fMatSP,rcotr(GaSP));
+ t1 = rcoatb(rcotr(fMatSP),rcotr(GaSP));
fMatSP = GaSP = [];
AtA = rcoata(GxSP);
@@ -3147,8 +3146,8 @@ func aoinit(disp=,clean=,forcemat=,svd=,dpi=,keepdmconfig=)
DtermSP = rcoadd(t1,GxSP);
t1 = [];
(*GxSP.xn) *= -1;
+ t2 = rcoatb(DtermSP,GxSP);
- t2 = rcotr(rcoatb(GxSP,DtermSP));
DtermSP = [];
CphiSPrco = ruo2rco(CphiSP);
CphiSP = [];
@@ -3812,7 +3811,7 @@ func go(nshot,all=)
(*loop.gainho)(order-1) * dm(nm).gain * errmb(indexDm(1,nm):indexDm(2,nm),imb);
}
} else { // tomographic DM; DM commands from virtual DMs
- virtualDMs = *dm(nm).fitvirtualdm;
+ virtualDMs = int(*dm(nm).fitvirtualdm);
virtualdmcommand = [];
for (idx=1;idx<= numberof(virtualDMs);idx++){
grow, virtualdmcommand, *dm(virtualDMs(idx))._command;
@@ -3832,7 +3831,7 @@ func go(nshot,all=)
yv = *dm(nm)._y - avg(*dm(nm)._y);
yv = yv / sqrt(sum(yv*yv));
-
+
*dm(nm)._command -= avg(*dm(nm)._command);
*dm(nm)._command -= (xv(+)*(*dm(nm)._command)(+))*xv;
*dm(nm)._command -= (yv(+)*(*dm(nm)._command)(+))*yv;
@@ -3904,7 +3903,7 @@ func go(nshot,all=)
okcscreen = ( is_set(controlscreen) && (((i-1) % controlscreen) == 0) );
if (is_set(controlscreen) && (i == loop.niter)) okcscreen=1; // display at last iteration
- if (savephase||okdisp) {
+ if (okdisp) {
// get the residual phase; initially for the first target
// display and save the residual wavefront
residual_phase=get_phase2d_from_dms(1,"target") +
@@ -3918,8 +3917,6 @@ func go(nshot,all=)
if (savephase){result=yao_fitswrite(YAO_SAVEPATH+"/"+parprefix+"_rwf"+swrite(loopCounter,format="%i")+".fits",float(residual_phase*pupil));}
}
-
-
// Computes the instantaneous PSF:
if (ok) {
if ((sim.svipc>>1)&1) { // use child

0 comments on commit 267aae8

Please sign in to comment.