~ubuntu-branches/ubuntu/trusty/yorick-spydr/trusty : revision 1

1

/*

2

* spydr_psffit.i

3

*

4

* PSF fitting functions for spydr. Computes Strehl ratio and FWHM.

5

*

6

* This file is part of spydr, an image viewer/data analysis tool

7

*

8

* $Id: spydr_psffit.i,v 1.10 2008/02/12 13:58:43 frigaut Exp $

9

*

10

11

*

12

* This program is free software; you can redistribute it and/or modify it

13

* under the terms of the GNU General Public License as published by the

14

* Free Software Foundation; either version 3 of the License, or (at your

15

* option) any later version.

16

*

17

* This program is distributed in the hope that it will be useful, but

18

* WITHOUT ANY WARRANTY; without even the implied warranty of

19

* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU

20

* General Public License for more details.

21

*

22

* You should have received a copy of the GNU General Public License

23

* along with this program. If not, see <http://www.gnu.org/licenses/>.

24

*

25

* $Log: spydr_psffit.i,v $

26

* Revision 1.10 2008/02/12 13:58:43 frigaut

27

* changelog to version 0.7.7:

28

*

29

* - fixed a bug when spydr_lut is not 0 and one creates a new

30

* window.

31

* - other minor bug fixes.

32

* - updated spydr man page

33

* - written and published web doc on maumae.

34

*

35

* Revision 1.9 2008/02/10 15:08:07 frigaut

36

* Version 0.7.6:

37

* - can now change the dpi on the fly. ctrl++ and ctrl+- will enlarge

38

* or shrink the graphical areas. long time missing in yorick.

39

* I have tried to make the window resizable, but it's a mess. Not

40

* only in the management of events, but also in the policy: really,

41

* only enlarging proportionally makes sense.

42

* - changed a bit the zoom behavior: now zoom is started once (the first

43

* time the mouse enter drawingarea1), and does not stop from that point.

44

* This is not ideal/economical (although disp_zoom returns immediately

45

* if the mouse is not in the image window), but it has the advantage

46

* of being sure the disp_zoom process does not spawn multiple instances

47

* (recurrent issue with "after").

48

* - The menu items in the left menu bar are hidden/shown according to the

49

* window size.

50

* - gotten rid of a few (unused) functions in spydr.i (the progressbar

51

* and message functions) that were conflicting with other pyk instances.

52

* - there's now focus in and out functions that will reset the current

53

* window to what it was before the focus was given to spydr. This is

54

* convenient when one just want to popup a spydr window to look at an

55

* image, and then come back to whatever one was doing without having to

56

* execute a window,n command.

57

* - fixed a bug in disp_cpc. Now, when a "e"/"E" command is executed

58

* while a subimage is displayed, the "e"/"E" applies to the displayed

59

* subimage, not the whole image.

60

* - changed a bit the behavior of the lower graphical area: not the y

61

* range is the same as the image zcuts (cmin/cmax).

62

* - fixed a small bug in get_subim (using floor/ceil instead of round

63

* for the indices determination).

64

* - added "compact" keyword to the spydr function (when called from

65

* within yorick).

66

* - clipping dpi values to [30,400].

67

* - spydr.py: went for a self autodic instead of an explicit

68

* declaration of all functions.

69

* - implemented smoothing by x2

70

* - implemented 1d linear fitting

71

*

72

* Revision 1.8 2008/02/02 05:12:05 frigaut

73

* fixed bug when picking star for fitting while being in "graphical axis

74

* in arcsec" mode.

75

*

76

* Revision 1.7 2008/01/30 05:28:19 frigaut

77

* - added spydr_pyk to avoid conflicts with other calls of pyk, and modify

78

* spydr_pyk for our purpose. I know this means we will not benefit from

79

* future pyk code improvements, but I can deal with that.

80

* - added check of yorick main version to avoid use with V<2.1.05 (in which

81

* current_mouse does not exist)

82

*

83

* Revision 1.6 2008/01/25 03:03:49 frigaut

84

* - updated license or license text to GPLv3 in all files

85

*

86

* Revision 1.5 2008/01/24 15:05:17 frigaut

87

* - added "delete from stack" feature

88

* - some bugfix in psffit

89

*

90

* Revision 1.4 2008/01/23 21:11:22 frigaut

91

* - load of new things:

92

*

93

* New Features:

94

* - added a number of command line flags (see man page or spydr -h)

95

* - can now handle series of image of different sizes

96

* - can mix single image and cube

97

* - cmin and cmax are now set per image (sticky setting)

98

* - image titles are better handled

99

* - updated man page

100

* - new image can be opened from the GUI menu (filechooser, multiple

101

* selection ok)

102

* - migrated to a spydrs structure, replaced many different variables, cleaner.

103

* - now opens the GUI even with no image argument (can use "open" from menu)

104

* - all errors are now also displayed as popups (critical quits yorick

105

* when called from shell)

106

* - because some (of the more critical) errors can happen before python is

107

* started, I had to use zenity for the popup window. New dependency.

108

* - added an "append" keyword to spydr. If set, the new image is appended

109

* to the list of displayed image. The old ones are kept, and the total

110

* number of image is ++

111

* - append is also available from the GUI menu

112

* - any action on displayed image can be null by using "help->refresh

113

* display" (in particular, sigmafilter)

114

* - created "about" dialog.

115

* - added an "image" menu (with names of all images in stack). user can

116

* select image form there.

117

* - added an "ops" (operation) menu. Can compute median, average, sum and

118

* rms of cube.

119

* - small gui (without lower panel) form is called with --compact (-c)

120

*

121

* Bug fixes:

122

* - fixed path to find python and glade files

123

* - fixed path for configuration file

124

* - main routine re-written and much more robust and clean

125

* - (kind of) solved a issue where image got displayed several times

126

* because of echo from setting cmin and cmax

127

* - fixed thibaut bug when closing window.

128

* - fixed "called_from_shell" when no image argument.

129

* - waiting for a doc for the user buttons, set to insivible.

130

* - waiting for a proper implementation of find, pane set to invisible.

131

*

132

*

133

* - bug: sometimes the next/previous image does not register

134

*

135

* Revision 1.3 2008/01/17 13:15:17 frigaut

136

* - modified the name of lmfit (-> spydr_lmfit) in spydr_psffit to avoid

137

* conflicts with the lmfit of yutils.

138

* - modified all calls of lmfit -> spydr_lmfit (spydr_psffit and spydr_various)

139

* - also default box size is now 51 (was 181)

140

*

141

* Revision 1.2 2007/12/13 13:43:27 frigaut

142

* - added license headers in all files

143

* - added LICENSE

144

* - slightly modified Makefile

145

* - updated info

146

* - bumped to 0.5.1

147

*

148

*

149

*

150

*/

151

152

153

version = "1.5";

154

modifDate = "June 17, 2007";

155

156

require,"random.i";

157

require,"string.i";

158

require,"random_et.i";

159

require,"utils.i";

160

require,"spydr_various.i";

161

require,"astro_util1.i"; // for sky()

162

require,"aoutil.i"; // for fwhmStrehl()

163

164

struct s_yfwhmres { double xpos, xposerr, ypos, yposerr, pstrehl, pfwhm, xfwhm, xfwhmerr, yfwhm, yfwhmerr, flux, fluxerr, el, elerr, angle, maxim, background;};

165

166

func parseDate(strdate)

167

{

168

y=m=d=0;

169

sread,strdate,format="%4d-%2d-%2d",y,m,d;

170

return [y,m,d];

171

}

172

func parseTime(strtime)

173

{

174

h=m=0; s=0.;

175

sread,strtime,format="%2d:%2d:%f",h,m,s;

176

return h+m/60.+s/3600.;

177

}

178

179

func getam(hdr,verbose=)

180

/* DOCUMENT func getam(hdr,verbose=)

181

Determine airmass from the information available in the AC header

182

for use in yfwhm.i procedure.

183

SEE ALSO:

280

*/

281

{

282

extern lmfit_amin,lmfit_amax;

283

w = where(lmfit_amin!=lmfit_amax);

284

if (numberof(w)!=0) {

285

if (dir==1) {

286

a(w) = (atan(a(w))/pi+0.5)*(lmfit_amax(w)-lmfit_amin(w))+lmfit_amin(w);

287

} else if (dir==-1) {

288

a(w) = tan(((a(w)-lmfit_amin(w))/(lmfit_amax(w)-lmfit_amin(w))-0.5)*pi);

289

} else error,"dir has to be 1 (forward) or -1 (backward) ";

290

}

291

return a;

292

}

293

294

func gaussian(x,ai)

295

{

296

local a; a=ai;

297

298

// a = [sky,Totalflux,Xcent,Ycent,~Xfwhm,~Yfwhm,angle]

299

a = lmfit_lim(a,1);

300

alpha = a(7)/180.*pi;

301

xp = (x(,,1)-a(3))*cos(alpha)+(x(,,2)-a(4))*sin(alpha);

302

yp = -(x(,,1)-a(3))*sin(alpha)+(x(,,2)-a(4))*cos(alpha);

303

r = sqrt(xp^2.+yp^2.);

304

if (a(5)==0) {

305

z = exp(-(yp/a(6))^2.);

306

} else if (a(6)==0) {

307

z = exp(-(xp/a(5))^2.);

308

} else z = exp(-((xp/a(5))^2.+(yp/a(6))^2.));

309

if (sum(z)==0) return a(1)+z;

310

z = a(1)+a(2)*z/sum(z);

311

return z;

312

}

313

314

315

func moffatRound(x,a)

316

{

317

// a=[sky,total,xc,yc,a,coefpow]

318

// a6 = atan(a(6))/(pi/2.)+1.2;

319

a1 = a(5);

320

a(6)=clip(a(6),-30,30);

321

xp = x(,,1)-a(3);

322

yp = x(,,2)-a(4);

323

// z = (1. + ((xp/a1)^2.+(yp/a1)^2.))^(-a(6));

324

if (a(6)==0) {

325

z = 1.+zp*0.;

326

} else {

327

if (a1==0) {

328

z = (1. + ((yp/a1)^2.))^(-a(6));

329

} else {

330

z = (1. + ((xp/a1)^2.+(yp/a1)^2.))^(-a(6));

331

}

332

}

333

if (sum(z)==0) return a(1)+z;

334

z = a(1)+a(2)*z/sum(z);

335

return z;

336

}

337

338

func moffat(x,ai)

339

{

340

local a; a=ai;

341

// a=[sky,total,xc,yc,a,b,angle,coefpow]

342

// alpha angle of longest axis clockwise

343

// a(8)=clip(a(8),-50,50);

344

a = lmfit_lim(a,1);

345

alpha = a(7)/180.*pi;

346

a1 = a(5);

347

a2 = a(6);

348

xp = (x(,,1)-a(3))*cos(alpha)+(x(,,2)-a(4))*sin(alpha);

349

yp = -(x(,,1)-a(3))*sin(alpha)+(x(,,2)-a(4))*cos(alpha);

350

if (a(8)==0) {

351

z = 1.+zp*0.;

352

} else {

353

if (a1==0) {

354

z = (1. + ((yp/a2)^2.))^(-a(8));

355

} else if (a2==0) {

356

z = (1. + ((xp/a1)^2.))^(-a(8));

357

} else {

358

z = (1. + ((xp/a1)^2.+(yp/a2)^2.))^(-a(8));

359

}

360

}

361

if (sum(z)==0) return a(1)+z;

362

z = a(1)+a(2)*z/sum(z);

363

return z;

364

}

365

366

func printhelp(void)

367

{

368

write,"Yorick function to interactively measure FWHM on an image";

369

write,

370

"Syntax: yfwhm [-help] [-a -p pixelsize -b boxsize -f functype -e extnum -s saturation -mag -v] image";

371

}

372

func printlonghelp(void)

373

{

374

write,"Yorick function to interactively measure FWHM on an image";

375

write,"Version "+version+" / Last modified "+modifDate;

376

write,"Syntax: yfwhm [-help] [-p pixelsize -b boxsize -mag] image";

377

write,"-help Prints this message";

378

// write,"-1 Uses only one window for graphic display";

379

write,"-a Outputs airmass corrected FWHM values";

380

write,"-p pixelsize Specify the image pixel size";

381

write,"-b boxsize Specify the size of the box of sub-images,";

382

write," usually 4-10 times the fwhm";

383

write,"-f functype function to use for fit (gaussian,special,moffat)";

384

write,"-e extnum fits extension number (main = 1)";

385

write,"-s saturation Saturation value (prevent picking saturated stars)";

386

write,"-mag Output flux in magnitude (zp=spydr_zero_point is used)";

387

write,"-v Verbose mode (more chatty)";

388

}

389

390

func yfwhm(bim,onepass,xstar,ystar,fluxstar,boxsize=,saturation=,pixsize=,funtype=, \

391

magswitch=,verbose=,airmass=)

392

/* DOCUMENT func yfwhm(image,boxsize=,saturation=,pixsize=,funtype=,

393

magswitch=,nwindow=,verbose=,airmass=)

394

image = 2D image

395

airmass = airmass. Outputs airmass corrected FWHM values

396

pixsize = Specify the image pixel size

397

boxsize = Specify the size of the box of sub-images

398

(usually 4-10 times the fwhm)

399

funtype = function to use for fit (gaussian,special,moffat

400

saturation = Saturation value (prevents picking saturated stars)

401

magswitch = Output flux in magnitude (zp=spydr_zero_point is used)

402

nwindow = Number of window for UI (default 2)

403

verbose = Verbose mode (0/1=more chatty)

404

x and y: added 2007jun15 for compatibilty with spydr find mode.

405

if x and y are set, then the interactive mode is turned off,

406

and the (x,y) coordinates are looped on to produce the

407

final yfwhmres (this function loops on the (x,y) and fit

408

each image in turn).

409

*/

410

{

411

extern spydr_fit_fwhm_estimate,imnum;

412

extern spydr_fit_background_estimate;

413

414

if (!is_set(boxsize)) boxsize = spydr_boxsize;

415

if (!is_set(saturation)) saturation = spydr_saturation;

416

if (!is_set(pixsize)) {pixsize = spydrs(imnum).pixsize;}

417

if (pixsize!=1.0) pixset=1;

418

if (!is_set(funtype)) {funtype = spydr_funtype;} else {funcset=1;};

419

if (!is_set(magswitch)) magswitch= output_magnitudes;

420

if (!is_set(airmass)) airmass = spydr_airmass;

421

nwindow=2;

422

if (!is_set(verbose)) verbose=0;

423

424

if (!compute_strehl) {

425

if (funtype == "gaussian") {write,"Using Gaussian fit";}

426

if (funtype == "special") {write,"Using Special fit";}

427

if (funtype == "moffat") {write,"Using Moffat fit";}

428

}

429

430

show_lower_gui,1;

431

432

batch_mode = (xstar!=[]);

433

434

if (batch_mode) {

435

if (numberof(xstar)!=numberof(ystar)) \

436

error,"X and Y are set but do not have the same dim";

437

n_to_do = numberof(xstar);

438

write,"Entering non interactive mode";

439

}

440

441

// spydr_disp;

442

window,spydr_wins(1);

443

444

maskarg = array(1,narg);

445

446

yfwhmres = s_yfwhmres();

447

allres = [];

448

449

local b,pow,zp,f,ferr,el,eler,an,airmass,dims,sky1,bim;

450

451

b = boxsize/2;

452

// update zoom to match boxsize:

453

rad4zoom = b;

454

pow = 0.85;

455

zp = spydr_zero_point;

456

f = array(float,2,1);

457

ferr = array(float,2,1);

458

el = 0.;

459

eler = 0.;

460

an = 0.;

461

airmass = double(airmass);

462

463

dims = (dimsof(bim))(2:3);

464

465

sky1 = sky(bim,dev1);

466

bim = bim-sky1;

467

if (saturation != 0.) {saturation -= sky1;}

468

469

/* this is called from spydr, so the windows pre-exist. */

470

471

if (onepass) {

472

// write,"Click on star";

473

} else {

474

write,"Left click on star for FWHM. Right click to exit.";

475

write,"Middle click to remove last entry.";

476

}

477

478

if (onepass) spydr_pyk_status_push,"Click on star";

479

else spydr_pyk_status_push,"BUTTONS: Left:Select Star / Middle:Remove last entry / Right:Exit.";

480

481

if (!compute_strehl) {

482

if (pixset) {

483

if (!magswitch) {

484

write,"X[pix] Y[pix] X FWHM[\"] Y FWHM[\"] FLUX[ADU] ELLIP ANGLE MAX";

485

} else {

486

write,"X[pix] Y[pix] X FWHM[\"] Y FWHM[\"] MAGNITUDE ELLIP ANGLE MAX";

487

}

488

} else {

489

if (!magswitch) {

490

write,"X[pix] Y[pix] X FWHM[pix] Y FWHM[pix] FLUX[ADU] ELLIP ANGLE MAX";

491

} else {

492

write,"X[pix] Y[pix] X FWHM[pix] Y FWHM[pix] MAGNITUDE ELLIP ANGLE MAX";

493

}

494

}

495

}

496

497

local nloop;

498

nloop=1;

499

500

do {

501

if (batch_mode) {

502

c = long(_(xstar(nloop),ystar(nloop)));

503

if (nloop==n_to_do) but=3; else but=1; // but=3 will exit main loop

504

} else { // interactive mode

505

res = mouse(1,0,"");

506

if (spydr_plot_in_arcsec) res(1:4) = spydr_arcsec_to_pixels(res(1:4));

507

spydr_pyk_status_push,"Processing...";

508

c = long(res(1:2));

509

but = res(10);

510

if (but == 3) break;

511

if (but == 2) {

512

if (numberof(el) == 1) {

513

write,"You can only unbuffer after having buffered at least one star!";

514

spydr_pyk_warning,"You can only unbuffer after having buffered at least one star!";

515

continue;

516

}

517

f = f(,:-1);

518

ferr = ferr(,:-1);

519

el = el(:-1);

520

eler = eler(:-1);

521

an = an(:-1);

522

write,"Last measurement taken out of star list";

523

spydr_pyk_warning,"Last measurement taken out of star list";

524

continue;

525

}

526

}

527

528

i1 = clip(c(1)-b,1,);

529

i2 = clip(c(1)+b,,dims(1));

530

j1 = clip(c(2)-b,1,);

531

j2 = clip(c(2)+b,,dims(2));

532

533

im = smooth(bim(i1:i2,j1:j2),2);

534

wm = where2(im == max(im))(*)(1:2)-b-1;

535

c = c + wm;

536

im = bim(i1:i2,j1:j2);

537

pos = c(1:2)-b;

538

pos = [i1,j1]-1;

539

// im = sigmaFilter(im,5,iter=2,silent=1);

540

if ((saturation > 0) && (max(im) > saturation)) {

541

if (onepass) {

542

write,"Some pixels > specified saturation level. Aborting !";

543

spydr_pyk_status_push,"Some pixels > specified saturation level. Aborting !";

544

exit;

545

} else {

546

write,"Some pixels > specified saturation level. Choose another star";

547

spydr_pyk_status_push,"Some pixels > specified saturation level. Choose another star";

548

continue;

549

}

550

}

551

sky2 = sky(im,dev2);

552

im = im - sky2;

553

d = dimsof(im);

554

555

w = 1.+0.*clip(im,dev2,)^2;

556

557

x = indices(d);

558

559

if (catch(0x11)) {

560

if (onepass) {

561

write,"Error detected, exiting.";

562

spydr_pyk_status_push,"Error detected, exiting.";

563

return;

564

}

565

write,"Error detected, skipping source";

566

spydr_pyk_status_push,"Error detected, skipping source";

567

nloop++;

568

continue;

569

}

570

571

extern lmfit_amin,lmfit_amax;

572

573

if (compute_strehl==0) {

574

if (funtype == "gaussian") {

575

// a = [sky,Totalflux,Xcent,Ycent,fwhm_parameter]

576

ai = [0,sum(im-median(im(*))),d(2)/2.,d(3)/2.,4.];

577

if (batch_mode) { // we've been given (guessed) coordinates, use them

578

ai = [0,fluxstar(nloop)*10.,xstar(nloop)-i1+1,ystar(nloop)-j1+1,spydr_fit_fwhm_estimate];

579

}

580

//r=spydr_lmfit(gaussianRound,x,ai,im,w,tol=1e-6,itmax=50,silent=1);

581

//ai(5)=abs(ai(5);

582

//a=[sky,Totalflux,Xcent,Ycent,a,b,angle]

583

a = [ai(1),ai(2),ai(3),ai(4),ai(5),ai(5),0.];

584

if (batch_mode) {

585

lmfit_amin = [0,0,a(3)-3,a(4)-3,a(5)-2,a(6)-2,0];

586

lmfit_amax = [0,0,a(3)+3,a(4)+3,a(5)+2,a(6)+2,0];

587

a(3:6) *=0.;

588

} else { // interative, no constraints.

589

lmfit_amin = [0,0,0,0,0,0,0];

590

lmfit_amax = [0,0,0,0,0,0,0];

591

}

592

r = spydr_lmfit(gaussian,x,a,im,w,stdev=1,tol=1e-8,itmax=50,silent=1);

593

tmp = gaussian(x,a);

594

a = lmfit_lim(a,1);

595

a(5:6) = abs(a(5:6));

596

err = *r.stdev;

597

pos = pos + a(3:4)-0.5;

598

angle = (a(7) % 180.) ;

599

angle -= 90.; // w.r.t. vertical axis instead of horizontal axis.

600

if (a(5) < a(6)) {

601

tmp2 = a(5); a(5) = a(6); a(6) = tmp2;

602

angle += 90;

603

}

604

while (angle > 90) angle -= 180.;

605

while (angle < -90) angle += 180.;

606

//if (angle < 0) {angle = angle+180.;}

607

//if (a(5) < a(6)) {angle = angle+90;}

608

//angle = (angle % 180.) ;

609

fwhm = a(5:6)*2*(-log(0.5))^(1./2.)*pixsize; //gaussian

610

fwhmerr = err(5:6)*2*(-log(0.5))^(1./2.)*pixsize;

611

fwhm = fwhm/airmass^0.6; fwhmerr = fwhmerr/airmass^0.6;

612

ellip = abs(fwhm(2)-fwhm(1))/avg(fwhm);

613

ellerr= 2*(fwhmerr(1)+fwhmerr(2))*(2*fwhm(2))/(fwhm(1)+fwhm(2))^2.;

614

615

} else if (funtype == "moffat") {

616

617

// a = [sky,total,xc,yc,fwhm_parameter,beta]

618

ai = [0,sum(im-median(im(*))),d(2)/2.,d(3)/2.,5.,1.7];

619

if (batch_mode) { // we've been given (guessed) coordinates, use them

620

if (spydr_fit_fwhm_estimate==[]) spydr_fit_fwhm_estimate=3.;

621

ai = [0,fluxstar(nloop)*10.,xstar(nloop)-i1+1,ystar(nloop)-j1+1,

622

clip(spydr_fit_fwhm_estimate,2.1,),1.7];

623

}

624

// r = spydr_lmfit(moffatRound,x,ai,im,w,tol=1e-6,itmax=50,silent=1);

625

// a = [sky, total, xc, yc, a, b, angle,beta]

626

a = [ai(1),ai(2),ai(3),ai(4),ai(5),ai(5),0.,ai(6)];

627

// a = [ai(1),ai(2),ai(3),ai(4),ai(5),ai(5),0.];

628

if (batch_mode) {

629

lmfit_amin = [0,0,a(3)-3,a(4)-3,a(5)-2,a(6)-2,0,1.];

630

lmfit_amax = [0,0,a(3)+3,a(4)+3,a(5)+2,a(6)+2,0,3.];

631

a(3:6) *=0.;

632

} else { // interative, no constraints.

633

lmfit_amin = [0,0,0,0,0,0,0,0];

634

lmfit_amax = [0,0,0,0,0,0,0,0];

635

}

636

fit = [1,2,3,4,5,6,7,8];

637

r = spydr_lmfit(moffat,x,a,im,w,stdev=1,tol=2e-8,itmax=50,silent=1,fit=fit);

638

tmp = moffat(x,a);

639

a = lmfit_lim(a,1);

640

a(5:6) = abs(a(5:6));

641

err = *r.stdev;

642

pos = pos + a(3:4)-0.5;

643

angle = (a(7) % 180.) ; //write,format="%f %f %f\n",a(5),a(6),angle;

644

angle -= 90.; // w.r.t. vertical axis instead of horizontal axis.

645

if (a(5) < a(6)) {

646

tmp2 = a(5); a(5) = a(6); a(6) = tmp2;

647

angle += 90;

648

}

649

while (angle > 90) angle -= 180.;

650

while (angle < -90) angle += 180.;

651

652

fwhm = 2*a(5:6)*sqrt(0.5^(-1./a(8))-1.)*pixsize; // moffat

653

fwhmerr = fwhm*(err(5:6)/a(5:6)+

654

0.5*abs(log(0.5))*err(8)/a(8)^2.*0.5^(1./a(8))/(0.5^(1./a(8))-1.));

655

fwhm = fwhm/airmass^0.6; fwhmerr = fwhmerr/airmass^0.6;

656

fwhmerr = fwhmerr(sort(fwhm)(::-1));

657

fwhm = fwhm(sort(fwhm)(::-1));

658

ellip = (fwhm(1)-fwhm(2))/avg(fwhm);

659

ellerr = 2*(fwhmerr(1)+fwhmerr(2))*(2*fwhm(2))/(fwhm(1)+fwhm(2))^2.;

660

}

661

662

maxim = max(tmp);

663

window,spydr_wins(3);

664

tv,transpose(grow(transpose(im),transpose(tmp),

665

transpose(im-tmp+a(1)))),square=1;

666

plt,"image",0.201,0.89,tosys=0;

667

plt,"fit",0.357,0.89,tosys=0;

668

plt,"image-fit",0.51,0.89,tosys=0;

669

spydr_xytitles,"pixels","pixels";

670

window,spydr_wins(1);

671

672

grow,f,fwhm;

673

grow,ferr,fwhmerr;

674

grow,el,ellip;

675

grow,eler,ellerr;

676

grow,an,angle;

677

678

if (magswitch) {flux = zp-2.5*log10(clip(a(2),1e-10,));} else {flux = a(2);}

679

680

yfwhmres.xpos = pos(1);

681

yfwhmres.ypos = pos(2);

682

yfwhmres.xposerr = (*r.stdev)(3);

683

yfwhmres.yposerr = (*r.stdev)(4);

684

yfwhmres.xfwhm = fwhm(1);

685

yfwhmres.yfwhm = fwhm(2);

686

yfwhmres.xfwhmerr = fwhmerr(1);

687

yfwhmres.yfwhmerr = fwhmerr(2);

688

yfwhmres.flux = flux;

689

yfwhmres.fluxerr = (*r.stdev)(2);

690

yfwhmres.el = ellip;

691

yfwhmres.elerr = ellerr;

692

yfwhmres.angle = angle;

693

yfwhmres.maxim = maxim;

694

yfwhmres.background = a(1)+sky1+sky2;

695

696

grow,allres,yfwhmres;

697

698

msg=swrite(format="%7.2f %7.2f %6.3f+/-%5.3f %6.3f+/-%5.3f %9.1f %4.2f %6.2f %6.1f",

699

pos(1),pos(2),fwhm(1),fwhmerr(1),fwhm(2),fwhmerr(2),flux,ellip,angle,maxim);

700

write,format="%s\n",msg;

701

702

703

pos(1),pos(2),fwhm(1),fwhm(2),flux,ellip,angle,maxim,a(1)+sky1+sky2);

704

msg = msg+" ("+funtype+")";

705

706

} else if (compute_strehl) {

707

708

if (pixset==0) {

709

spydr_pyk_error,"Need pixsize set to compute Strehl!";

710

return;

711

}

712

if (spydrs(imnum).wavelength==0) {

713

spydr_pyk_warning,"Need wavelength to compute Strehl!";

714

return;

715

}

716

717

// determination of background

718

sdim = dimsof(im)(2);

719

rmask = spydr_strehlaper/2.;

720

// write,format="rmask=%f, boxsize=%f\n",rmask*1.,boxsize*1.;

721

if ((2*rmask)>boxsize) {

722

spydr_pyk_error,"boxsize smaller than aperture for Strehl. Increase boxsize or decrease aperture.";

723

return;

724

}

725

smask = (dist(sdim)>rmask);

726

psky = im(where(smask)); // outside of disk

727

nsig=4.;

728

// first pass:

729

w = where( (psky>(median(psky)-psky(rms)*nsig)) & (psky<(median(psky)+psky(rms)*nsig)) );

730

skyavg = psky(w)(avg);

731

skyrms = psky(w)(rms);

732

// second pass:

733

w = where( (psky>(skyavg-skyrms*nsig)) & (psky<(skyavg+skyrms*nsig)) );

734

hy = histo2(psky(w),hx);

735

ag = [max(hy),hx(wheremax(hy)(1)),3.];

736

clmfit,float(hy),hx,ag,"a(1)*exp(-0.5*((x-a(2))/a(3))^2.)",yfit;

737

skyavg = ag(2);

738

skyrms = ag(3)*2.355;

739

// plots:

740

curw = current_window();

741

window,spydr_wins(3);

742

fma; limits,square=0; limits;

743

plh,hy,hx;

744

plh,yfit,hx,color="red";

745

spydr_pltitle,swrite(format="Background and fit (avg=%f, rms=%f)",skyavg,skyrms);

746

spydr_xytitles,"value","number in bin";

747

window,curw;

748

//write,format="a(1) = %f ; sky = %f +/- %f\n",a(1),skyavg,skyrms;

749

750

fwhmStrehl,im-skyavg,pixsize,spydrs(imnum).wavelength,spydr_teldiam, \

751

spydr_cobs,pstrehl,pfwhm,rmask=rmask,silent=1,source=spydr_sourcediam;

752

753

// apply fudge

754

pstrehl *= spydr_strehlfudge;

755

756

// plot boxes

757

758

if (spydr_plot_in_arcsec) {

759

fact = spydrs(imnum).pixsize;

760

axtit = "arcsec";

761

} else {

762

fact = 1.0f;

763

axtit = "pixels";

764

}

765

766

plg,([j2+1,j2+1,j1,j1,j2+1])*fact,([i1,i2+1,i2+1,i1,i1])*fact,color="red";

767

plt,"Sky",(i1)*fact,(j2+1)*fact,justify="LT",tosys=1,color="red",opaque=0;

768

tmp = span(0.,2*pi,100);

769

plg,((j2+j1)/2.+0.5+rmask*sin(tmp))*fact,((i2+i1)/2.+0.5+rmask*cos(tmp))*fact,color="red";

770

// print results

771

nptvalid = pi*rmask^2.;

772

nptsky = numberof(w);

773

strehl_err = skyrms/sqrt(nptsky)*nptvalid/max(im)*pstrehl;

774

// note: FIXME. This is just the strehl error due to the estimation

775

// of the average sky. There is an additional component

776

// due to the noise within the aperture.

777

msg = swrite(format="FWHM = %.3f | strehl = %.2f",pfwhm,pstrehl);

778

write,format="%s wvl=%.3f FWHM=%.3f Strehl=%.2f +/- %.2f (fudge=%.3f)\n",\

779

spydrs(imnum).name,spydrs(imnum).wavelength,pfwhm,pstrehl, \

780

strehl_err,spydr_strehlfudge;

781

//yfwhmres.pstrehl = pstrehl;

782

//yfwhmres.pfwhm = pfwhm;

783

}

784

785

spydr_pyk_status_push,msg;

786

787

if (onepass) break;

788

// typeReturn;

789

nloop++;

790

} while (but != 3);

791

792

793

if (!compute_strehl) {

794

f = f(,2:);

795

ferr = ferr(,2:);

796

el = el(2:);

797

eler = eler(2:);

798

avgfwhm = sum((f*1./ferr)(*))/sum(1./ferr(*));

799

// stdfwhm = f(*)(rms);

800

stdfwhm = avg([f(1,)(rms),f(2,)(rms)]); // avg X and Y rms

801

avgel = avg(el);

802

stdel = el(rms)+sqrt(sum(eler^2.))/numberof(eler);

803

804

if (pixset) {

805

msg=swrite(format="Median FWHM : X = %5.3f / Y = %5.3f / <XY> = %6.3f [arcsec]",

806

median(f(1,)),median(f(2,)),avg([median(f(1,)),median(f(2,))]));

807

} else {

808

msg=swrite(format="Median FWHM : X = %6.3f / Y = %6.3f / <XY> = %6.3f [pixel]",

809

median(f(1,)),median(f(2,)),avg([median(f(1,)),median(f(2,))]));

810

}

811

write,format="\n%s\n",msg;

812

if (!onepass) spydr_pyk_status_push,msg;

813

// in pixels:

814

spydr_fit_fwhm_estimate = avg([median(f(1,)),median(f(2,))])/pixsize;

815

spydr_fit_background_estimate = avg(allres.background);

816

} else { // compute_strehl

817

spydr_fit_fwhm_estimate = pfwhm/pixsize;

818

spydr_fit_background_estimate = skyavg;

819

}

820

821

spydr_pyk,swrite(format="y_text_parm_update('find_fwhm','%.3f')",spydr_fit_fwhm_estimate);

822

823

if (compute_strehl) return _(pfwhm,pstrehl);

824

else return allres;

825

}

826

827

828

require, "random.i";

829

830

struct lmfit_result {

831

/* DOCUMENT lmfit_result -- structure returned by lmfit

832

*/

833

long neval;

834

long niter;

835

long nfit;

836

long nfree;

837

long monte_carlo;

838

double chi2_first;

839

double chi2_last;

840

double conv;

841

double sigma;

842

double lambda;

843

pointer stdev;

844

pointer stdev_monte_carlo;

845

pointer correl;

846

};

847

848

func spydr_lmfit(f, x, &a, y, w, fit=, correl=, stdev=, gain=, tol=, deriv=, itmax=,

849

lambda=, eps=, monte_carlo=,silent=)

850

/* DOCUMENT spydr_lmfit

851

Non-linear least-squares fit by Levenberg-Marquardt method.

852

This is a local copy, slightly modified, of lmfit.i (in yutils).

853

For help, please refer to the lmfit document section.

854

In general, use lmfit instead of spydr_lmfit.

855

*/

856

{

857

local grad;

858

859

/* Maybe subset of parameters to fit. */

860

if (structof(a)!=double) {

861

a+= 0.0;

862

if (structof(a)!=double)

863

error, "bad data type for parameters (complex unsupported)";

864

}

865

na= numberof(a);

866

if (is_void(fit))

867

fit= indgen(na);

868

else if (dimsof(fit)(1) == 0)

869

fit= [fit];

870

nfit= numberof(fit);

871

if (!nfit)

872

error, "no parameters to fit.";

873

874

/* Check weights. */

875

if (is_void(w)) w= 1.0;

876

else if (anyof(w < 0.0))

877

error, "bad weights.";

878

if (numberof(w) != numberof(y))

879

w += array(0.0, dimsof(y));

880

nfree= sum(w != 0.0) - nfit; // Degrees of freedom

881

if (nfree <= 0)

882

error, "not enough data points.";

883

884

/* Other settings. */

885

diag= indgen(1:nfit^2:nfit+1); // Subscripts of diagonal elements

886

if (is_void(lambda)) lambda= 1e-3;

887

if (is_void(gain)) gain= 10.0;

888

if (is_void(itmax)) itmax= 100;

889

if (is_void(eps)) eps= 1e-6; // sqrt(machine_precision)/100

890

if (1.0+eps <= 1.0)

891

error, "bad value for EPS.";

892

if (is_void(tol)) tol= 1e-7;

893

monte_carlo= is_void(monte_carlo) ? 0 : long(monte_carlo);

894

warn_zero= 0;

895

warn= "*** Warning: spydr_lmfit ";

896

neval= 0;

897

conv= 0.0;

898

niter= 0;

899

900

while (1) {

901

if (deriv) {

902

m= f(x, a, grad, deriv=1);

903

neval++;

904

grad= nfit == na ? grad(*,) : grad(*,fit);

905

} else {

906

if (!niter) {

907

m= f(x, a);

908

neval++;

909

}

910

inc= eps * abs(a(fit));

911

if (numberof((i= where(inc <= 0.0)))) inc(i)= eps;

912

grad= array(double, numberof(y), nfit);

913

for (i=1; i<=nfit; i++) {

914

anew= a; // Copy current parameters

915

anew(fit(i)) += inc(i);

916

grad(,i)= (f(x,anew)-m)(*)/inc(i);

917

}

918

neval += nfit;

919

}

920

beta= w * (chi2= y-m);

921

if (niter) chi2= chi2new;

922

else chi2= chi2_first= sum(beta * chi2);

923

beta= grad(+,) * beta(*)(+);

924

alpha= ((w(*)(,-) * grad)(+,) * grad(+,));

925

gamma= sqrt(alpha(diag));

926

if (anyof(gamma <= 0.0)) {

927

/* Some derivatives are null (certainly because of rounding

928

* errors). */

929

if (!warn_zero) {

930

if (!silent) write, warn+"founds zero derivatives.";

931

warn_zero= 1;

932

}

933

gamma(where(gamma <= 0.0))= eps * max(gamma);

934

if (allof(gamma<=0.0)) goto done; // case where all gamma=0

935

}

936

gamma= 1.0 / gamma;

937

beta *= gamma;

938

alpha *= gamma(,-) * gamma(-,);

939

940

while (1) {

941

alpha(diag)= 1.0 + lambda;

942

anew= a;

943

anew(fit) += gamma * LUsolve(alpha, beta);

944

m= f(x, anew);

945

neval++;

946

d= y-m;

947

chi2new= sum(w*d*d);

948

if (chi2new < chi2)

949

break;

950

lambda *= gain;

951

if (allof(anew == a)) {

952

/* No change in parameters. */

953

if (!silent) write, warn+"makes no progress.";

954

goto done;

955

}

956

}

957

a= anew;

958

lambda /= gain;

959

niter++;

960

conv= 2.0*(chi2-chi2new)/(chi2+chi2new);

961

if (conv <= tol)

962

break;

963

if (niter >= itmax) {

964

if (!silent) {

965

write, format=warn+"reached maximum number of iterations (%d).\n",

966

itmax;

967

}

968

break;

969

}

970

}

971

972

done:

973

sigma= sqrt(nfree/chi2);

974

result= lmfit_result(neval=neval, niter=niter, nfree=nfree, nfit=nfit,

975

lambda=lambda, chi2_first=chi2_first, chi2_last=chi2, conv=conv,

976

sigma=sigma);

977

if (correl || stdev) {

978

/* Compute correlation matrice and/or standard deviation vector. */

979

alpha(diag)= 1.0;

980

alpha= LUsolve(alpha);

981

if (anyof((tmp1= alpha(diag)) < 0.0))

982

write, format=warn+"%s\n", "found negative variance(s)";

983

tmp1= sqrt(abs(tmp1));

984

if (stdev) {

985

/* Standard deviation is rescaled assuming that statistically

986

* chi2 = nfree +/- sqrt(2*nfree). */

987

(tmp2= array(double,na))(fit)= gamma * tmp1 / sigma;

988

result.stdev= &tmp2;

989

}

990

if (correl) {

991

gamma= 1.0 / tmp1;

992

alpha *= gamma(-,) * gamma(,-);

993

if (nfit == na) {

994

result.correl= α

995

} else {

996

(tmp2= array(double, na, na))(fit,fit)= alpha;

997

result.correl= &tmp2;

998

}

999

}

1000

}

1001

alpha= beta= gamma= []; // Free some memory.

1002

if (monte_carlo >= 1) {

1003

saa= 0.0*a;

1004

sig= (w > 0.0) /(sqrt(max(nfree/chi2*w, 0.0)) + (w == 0.0));

1005

for (i=1; i<=monte_carlo; i++) {

1006

anew= a;

1007

// ynew= y + sig * random_normal(dimsof(y));

1008

ynew= y + sig * random_n(dimsof(y));

1009

spydr_lmfit, f, x, anew, ynew, w, fit=fit, gain=gain, tol=tol,

1010

deriv=deriv, itmax=itmax, lambda=lambda, eps=eps;

1011

anew -= a;

1012

saa += anew * anew;

1013

}

1014

result.monte_carlo= monte_carlo;

1015

result.stdev_monte_carlo= &sqrt(saa / monte_carlo);

1016

}

1017

return result;

1018

}