1 scholl 1.4 c------------------------------------------------------------------------------------------------
2 c #define CODE_NAME "limbfit"
|
3 scholl 1.10 c #define CODE_VERSION "V5.2r0"
4 c #define CODE_DATE "Tue Mar 26 16:30:08 HST 2013"
|
5 scholl 1.4 c------------------------------------------------------------------------------------------------
|
6 arta 1.1 c Revision 1.0 2009/01/08 17:20:00 Marcelo Emilio
7 c changed assumed-size array declaration from "real xxx(1)" to "real xxx(3)"
8 c jreg=1024, jang=129,jprf=160, jpt=4000000, ahi=30000.0
9 c variable funcs changed to funcsb
10 c Revision 2.0 2010/18/12 11:00:00 Marcelo Emilio
11 c 2010/01/05 JK:Here's a version of limb_5.f that shouldn't return the center guess if clmt is
12 c input too large...lets try this at clmt of e-6, e-7, and e-8... -- J
13 c modified with input shift function to iteratively generate LDF
14
15 c Subroutine Limb
16 c Input:
17 c 1) Position and intensity values for limb annulus pixels in the real array
18 c "anls(3,*)" x,y coordinates: anls(1,*),anls(2,*); intensity anls(3,*).
19 c Expect the number of annulus pixels, integer "ndat", to be <= 40000
20 c 2) Intial guessimate of the center position, real "cmx,cmy"; try (511.5,
21 c 511.5)
|
22 scholl 1.8 c 3) Number of angular bins in which to find average profiles, integer "nang"
23 c to be <=16 (16 gives good resolution for a 7 pixel annulus)
|
24 arta 1.1 c 4) Number of radial bins, integer "nprf" <=40 (if this parameter is set to
|
25 scholl 1.8 c much less than 40, the order of the Chebyschev polynomial "jc" should be
26 c changed)
|
27 arta 1.1 c 5) Number of angular bins in which to compute alpha & beta, the limb scale
28 c factor and offset, integer "nreg" <= 512
29 c Output:
30 c 1) New improved center position
31 c 2) Mean radius of all pixels, real "rmean"
32 c 3) Center-finder convergence diagnostic: number iterations, integer "nitr"
33 c 4) LDF profile-fitting diagnostic: number of points cut, integer "ncut"
34 c 5) Average profile for "nang" angular bins: real array "lprf(nprf,nang+1)"
35 c where the last (nang+1) profile is the overall mean limb profile to which
36 c the Chebyschev polynomial is fit. The radial information is stored
37 c in real array "rprf(nprf)".
38 c 6) LDF scale factor and offset per angular bin, real arrays "alph(nreg)"
39 c and "beta(nreg)"
40 c 7) A general success/failure flag: integer "ifail". For successful
41 c operation ifail=0; any other value corresponds to a failure at a specific
42 c numbered site in the code which can be determined by searching for
43 c "ifail=#". With one exception, a failure (ifail >< 0) is coupled with
44 c a "return" to the limb wrapper calling procedure.
45 c 8) b0 is the "preshift" add this to the output beta to get the net
46 c *) centyp to specify if cmx/cmy are guess center(=0) or to be used as center(=1)
|
47 scholl 1.8 Subroutine limb(anls,npts,cmx,cmy,rguess,nitr,ncut,rprf,lprf,
48 & rsi, rso, dx, dy,
|
49 scholl 1.3 & alph,beta,ifail,b0, centyp, ahi)
|
50 arta 1.1
51 c parameter(jpt=40000,jreg=4096,jb=50,jc=50,jang=129,jprf=160)
|
52 scholl 1.6 c values for HMI rolls
53 implicit none
54 integer jpt,jb,jc,jreg,jang,jprf,nang,nprf,nreg
|
55 scholl 1.9 parameter(jpt=8000000,jb=200,jc=200,jreg=256,jang=181,jprf=64)
56 parameter(nang=180,nprf=64,nreg=256)
|
57 arta 1.1 integer i,j,itr,itr2,ilt,ilt2,lmt,ir,lista(3),setn,nremv,centyp
|
58 scholl 1.6 integer n,nb,ind,nc,spflag,nbad,wbad,cut,jnd
59 integer npts,nitr,ncut,ifail,nrem
60 real rguess,ain,aout,chebev
|
61 arta 1.1 real anls(3,jpt),an(3,jpt),cmx,cmy,alph(jreg),beta(jreg),rsi,rso
62 real pi,tpi,alo,ahi,flag,asi,aso,rni,rno,clmt,rem,dreg,dx,dy,fbrem
63 real rmax,rmin,rminp,rmaxp,dr,pix(jb+2),rad(jb+2),bin(jb+2)
64 real a,x,y,r,savr(jpt),savi(jpt),rmean,imean,rout,rin,theta
65 real inte(jpt),thet(jpt),sig(jpt),scmx,scmy,scrit,stpscl
|
66 scholl 1.6 real cf(3),cov(3,3),chi,cfo(3),cxo,cyo
67 real c(jc),cp(jc),cpp(jc),lsq,ln,dev,err,vr
|
68 arta 1.1 real d0,d1,d2,expn(6,jreg),da,db,dc,dd,de,df,dbb,daa
69 real crit,svx,svy,xscal,yscal
70 real lprf(jprf,jang),prf(jprf,jang),rprf(jprf),dprf,dang
71 real b0(jreg)
72 common /ldfcom/ spflag,nb,rad,bin
73 common /model/ nc,c,cp,cpp
74
|
75 scholl 1.6 external funcs
76 c print*,("in fortran")
77 c print*,("cc:"),centyp,ahi,rsi,rso,dx,dy
78 c print*,("--:"),ifail,ncut,nitr,cmx,cmy,rguess
79 c print*,("--:"),jprf,jang,jreg,jpt,npts
80 c print*,(">>"),rprf(1),rprf(jprf),lprf(1,1),lprf(jprf,jang)
81 c print*,(">>"),anls(1,1),anls(2,1),anls(3,1),anls(1,npts),anls(2,npts),anls(3,npts)
82 c print*,(">>:"),alph(1),alph(jreg),beta(1),beta(jreg),b0(1),b0(jreg)
|
83 arta 1.1 pi=3.141592654
84 tpi=2.0*pi
85 dreg=tpi/float(nreg) ! angular bin size for alpha & beta
86 dang=tpi/float(nang) ! angular bin size for profiles
|
87 scholl 1.6 alo=1.0 ! lower limit for acceptable data
|
88 scholl 1.3 c ahi=70000.0 ! upper limit ...
|
89 arta 1.1 flag=-2147483648.0 ! flag bad pixels
90 c rsi=18.0 ! rmean-rsi lower limit to use in center-finder
91 c rso=18.0 ! rmean+rso upper limit ...
92 asi=4.0 ! Imean*asi upper limit to use in center-finder
93 aso=0.5 ! Imean*aso lower limit ...
94 fbrem=2.0 ! number of std.dev. to cut data in sin/cos fit
95 c dx=-4.0 ! x increment to step by in center-finder
96 c dy=-4.0 ! y increment ...
97 stpscl=0.0001 ! limit on how small to scale the steps
98 ilt=10 ! "itr" loop cutoff for stepping dx,dy to center
99 ilt2=1 ! "itr2" cutoff for finding <r> w/ new center
100 c clmt=0.0000001 ! coefficient convergence criteria for center-finder
101 clmt=1E-7 ! coefficient convergence criteria for center-finder
102 rni=rsi ! as "rsi", but for remaining calculations
103 rno=rso ! as "rso", ...
104 nb=nprf ! number of radial bins (jb=50)
105 nc=16 ! order of the Chebyschev polynomial (jc=50)
106 rem=20.0 ! number of std.dev. to cut data in the Chebyschev fit
107 lmt=20 ! limit on # of cuts of rem*sigma; loop counter "cut"
108 c Initialize stuff; find the min & max radii
109 ifail=0
110 arta 1.1 svx=cmx
111 svy=cmy
112 do i=1,3
113 lista(i)=i
114 enddo
115 do i=1,jpt
116 inte(i)=0.0
117 thet(i)=0.0
118 savr(i)=0.0
119 savi(i)=0.0
120 enddo
121 rmax=0.0
122 rmin=100000.0
123 do 2 i=1,npts
124 a=anls(3,i)
125 if((a.le.alo).or.(a.gt.ahi)) goto 2
126 r=((anls(1,i)-cmx)**2+(anls(2,i)-cmy)**2)**0.5
127 if(r.gt.rmax) rmax=r
128 if(r.lt.rmin) rmin=r
129 ir=ir+1
130 2 enddo
131 arta 1.1 if(centyp.eq.1) goto 40
|
132 scholl 1.5 c PRINT*, ("center determination"),cmx,cmy
|
133 arta 1.1 c Begin center-finding routine. Starting from the rough center input, find
134 c the mean radius of the annulus. Define a sub-annulus; find intensity(angle).
135 c Iteratively, find a new center that minimizes the difference in the measured
136 c intensity and the function cf(1)*sin(theta)+cf(2)*cos(theta)+cf(3).
137 c The convergence criteria is cf(1)**2+cf(2)**2 <<< cf(3)**2. The code loops
138 c until this criteria is met or the counter "itr" reaches limit "ilt".
139 c After convergence, a new annulus and I(theta) is found with the new center.
140 c For robustness, this new array (with new mean radius) is tested for
141 c convergence until counter "itr2" reaches limit "ilt2"
142
143 itr2=0
144 scrit=1.0
145 nitr=0
146 10 continue ! loop ilt2 times s.t. test center w/ new radius
147 itr2=itr2+1
148
149 c Compute radial and intensity means
150 ir=0
151 do 20 i=1,npts
152 a=anls(3,i)
153 if((a.le.alo).or.(a.gt.ahi)) goto 20
154 arta 1.1 r=((anls(1,i)-cmx)**2+(anls(2,i)-cmy)**2)**0.5
155 ir=ir+1
156 savr(ir)=r
157 savi(ir)=a
|
158 scholl 1.6 20 enddo
159 c print*,ir
|
160 arta 1.1 call Hpsort(ir,savr)
161 call Hpsort(ir,savi)
162 i=ir/2
163 rmean=savr(i)
164 imean=savi(i)
165 c PRINT*, ("rmean rguess"), rmean, rguess
166
167 c Define annulus radial range
168 rin=rguess-rsi
169 rout=rguess+rso
170 c Define annulus with intensity range
171 ain=imean*asi
172 aout=imean*aso
173
174 itr=0
175 30 continue ! loop ilt times s.t. converge on a center
176 itr=itr+1
177
178 c PRINT*,"rmean, imean, ain, aout",rmean, imean, ain, aout
179 do i=1,jpt
180 sig(i)=1.0
181 arta 1.1 enddo
182 c Put sub-annulus data into intensity and angle arrays
183 call darr(anls,npts,ain,aout,cmx,cmy,rin,rout,thet,inte,setn)
184 c call gplot(setn, thet, inte, 1, hc, it)
185
186 c Fit sin/cos function (found in "funcs") to intensity(angle)
187 call lfit(thet,inte,sig,setn,cf,3,lista,3,cov,3,chi,funcs,ifail)
|
188 scholl 1.4 if(ifail.gt.0)then
189 c add by IS Oct5
190 ifail=20
191 return
192 endif
|
193 arta 1.1 c Remove outliers from fit
194 call fitb(thet,inte,setn,fbrem,chi,sig,cf,nremv)
195 c Fit again
196 call lfit(thet,inte,sig,setn,cf,3,lista,3,cov,3,chi,funcs,ifail)
|
197 scholl 1.4 if(ifail.gt.0)then
198 c add by IS Oct5
199 ifail=21
200 return
201 endif
|
202 arta 1.1 crit=(cf(1)*cf(1)+cf(2)*cf(2))/cf(3)/cf(3)
203 c Save the center with the best convergence
204 c PRINT*,"scrit,criti, cmx, cmy:",scrit,crit, cmx, cmy
205 if(crit.lt.scrit)then
206 scrit=crit
207 scmx=cmx
208 scmy=cmy
209 endif
210 c Test convergence; Save number of interations
211 if(crit.lt.clmt)then
212 nitr=nitr+itr
213 if(itr2.le.ilt2) then
214 c save file to debug
215 c OPEN(3,FILE='output.txt')
216 c WRITE(3,*)setn
217 c WRITE(3,'(I8, F18.8, F18.8)'),(j, thet(J), inte(j), J=1,setn)
218 c WRITE(3,*),(j, thet(J), inte(j),char(13)//char(10), J=1,setn)
219 c WRITE(3,*),(j, thet(J), inte(j), J=1,setn)
220
221 goto 10 ! loop again
222 endif
223 arta 1.1 c-1/2010 JRK added if condition to force one center calculation with big clmt
224 if(nitr.gt.2)goto 40 ! continue on only if one center correction
225 endif
226
227 c Recompute the coefficients for a center position offset by dx,dy
228 cfo(1)=cf(1)
229 cfo(2)=cf(2)
230 cxo=cmx
231 cmx=cmx+dx
232 cyo=cmy
233 cmy=cmy+dy
234 do i=1,jpt
235 sig(i)=1.0
236 enddo
237 call darr(anls,npts,ain,aout,cmx,cmy,rin,rout,thet,inte,setn)
238 call lfit(thet,inte,sig,setn,cf,3,lista,3,cov,3,chi,funcs,ifail)
|
239 scholl 1.4 if(ifail.gt.0)then
240 c add by IS Oct5
241 ifail=22
242 return
243 endif
|
244 arta 1.1 call fitb(thet,inte,setn,fbrem,chi,sig,cf,nremv)
245 call lfit(thet,inte,sig,setn,cf,3,lista,3,cov,3,chi,funcs,ifail)
246 c PRINT*, "ITR, CHI, NREMV", itr, chi, nremv
|
247 scholl 1.4 if(ifail.gt.0)then
248 c add by IS Oct5
249 ifail=23
250 return
251 endif
|
252 arta 1.1
253 c Compute the factor with two sets of coefficients to scale the step size
254 xscal=cfo(1)/(cf(1)-cfo(1))
255 yscal=cfo(2)/(cf(2)-cfo(2))
256 xscal=dx*xscal
257 yscal=dy*yscal
258 if ((abs(xscal).lt.stpscl).and.(abs(yscal).lt.stpscl)) goto 35
259 cmx=cxo-xscal
260 cmy=cyo-yscal
261 c PRINT*, "cmx, cmy, xscal, yscal, nitr", cmx, cmy, xscal, yscal, nitr
262 if(itr.le.ilt) goto 30 ! try new center
263 35 continue ! no convergence
264 if(itr2.le.ilt2) then
265 nitr=nitr+itr
266 goto 10 ! try new radius
267 endif
268 c No convergence; Reset center to that with the best criteria
269 c This is flagged in the output file by the parameter "nitr" being its maximum
270 c value and by "ifail"=1. Good data converges with under 5 iterations.
271 cmx=scmx
272 cmy=scmy
273 arta 1.1 nitr=nitr+itr
274 c Do not return here because convergence criteria could be too strict unless
275 c the center does not change or it changes alot
276 cxo=abs(cmx-svx)
277 cyo=abs(cmy-svy)
278 if(((cxo.eq.0.0).and.(cyo.eq.0.0)).or.
279 & ((cxo.gt.10.0).and.(cyo.gt.10.0)))then
280 ifail=1
281 return
282 endif
283
284 c Begin limb figure determination. Define a sub-annulus about the mean radius
285 c with the new center coordinates. Radially bin intensity data. Fit a Cheby-
286 c schev polynomial of order "nc". Iteratively, refit polynomial to data after
287 c outliers (> "rem"*sigma) have been removed. The limit to iteration counter
288 c "cut" is "lmt". Note: if cut > lmt, the code continues without action.
289
290 40 continue
291 c Find the min and max radii
292 rmax=0.0
293 rmin=100000.0
|
294 scholl 1.8 do 42 i=1,npts
295 a=anls(3,i)
296 if((a.le.alo).or.(a.gt.ahi)) goto 42
297 r=((anls(1,i)-cmx)**2+(anls(2,i)-cmy)**2)**0.5
298 if(r.gt.rmax) rmax=r
299 if(r.lt.rmin) rmin=r
300 42 enddo
|
301 arta 1.1
302 c Find a sub-annulus about the mean radius. This may not be necessary with the
303 c SOI data, but was for data where larger annuli (or full disk was available).
304 n=0
305 rin=rguess-rni
306 rout=rguess+rno
307 do 50 i=1,npts
308 a=anls(3,i)
309 if((a.le.alo).or.(a.gt.ahi)) goto 50
310 x=anls(1,i)
311 y=anls(2,i)
312 r=((x-cmx)**2+(y-cmy)**2)**0.5
313 if((r.lt.rin).or.(r.gt.rout)) goto 50
314 n=n+1
315 savr(n)=r
316 savi(n)=a
317 an(1,n)=x
318 an(2,n)=y
319 an(3,n)=a
320 50 enddo
321 c Loop to cut points deviating from the Chebyschev fit.
322 arta 1.1 nbad=0 ! current number of bad points
323 cut=0 ! loop counter
324 ncut=0 ! total points cut
325 60 continue
326
327 c Find min and max radii
328 rmax=0.0
329 rmin=100000.0
330 do 70 i=1,n
331 a=an(3,i)
332 if((a.le.alo).or.(a.gt.ahi)) goto 70
333 r=((an(1,i)-cmx)**2+(an(2,i)-cmy)**2)**0.5
334 if(r.gt.rmax) rmax=r
335 if(r.lt.rmin) rmin=r
|
336 scholl 1.6 70 enddo
|
337 arta 1.1 c Put data into "nb" radial bins
338 dr=(rmax-rmin)/float(nb)
339 rad(1)=rmin-dr/2.0
340 do i=2,nb+2
341 rad(i)=rad(i-1)+dr
342 c (IS: 1st bin -> rmin-dr , 2nd -> rmin)
343 enddo
|
344 scholl 1.6 c print*,rad
|
345 arta 1.1 do i=1,nb+2
346 pix(i)=0.0
347 bin(i)=0.0
348 enddo
349 nrem=0
|
350 scholl 1.6 c- print*,"-",alo,ahi,tpi,dreg,cmx,cmy,nrem
|
351 arta 1.1 do 80 i=1,n
352 a=an(3,i)
353 if((a.le.alo).or.(a.gt.ahi)) goto 80
354 nrem=nrem+1
355 r=((an(1,i)-cmx)**2+(an(2,i)-cmy)**2)**0.5
|
356 scholl 1.10 c this added to preshift limb distortion to get cleaner mean LDF
357 c (IS: maximum distorsion = 2 pixels)
|
358 arta 1.1 theta=atan2((an(2,i)-cmy),(an(1,i)-cmx))
359 if(theta.lt.0.0) theta=theta+tpi
360 ind=int(theta/dreg+1.0)
|
361 scholl 1.6 c- print*,"->",i,a,an(1,i),an(2,i),r,rmin,dr,ind,theta,b0(ind)
|
362 arta 1.1 r = r - b0(ind)
363 c end of radius correction
364 ind=int((r-rmin)/dr+2.0)
|
365 scholl 1.6 c- print*," ",r,rmin,dr,ind
|
366 scholl 1.4 c if(ind.lt.1) goto 80
367 c changed w/Jeff Oct 6,2011
368 if(ind.lt.1.or.ind.gt.nb) goto 80
|
369 arta 1.1 bin(ind)=bin(ind)+a
370 pix(ind)=pix(ind)+1.0
|
371 scholl 1.6 c print*,"=>",i,ind,bin(ind),pix(ind)
372 80 enddo
|
373 arta 1.1 do i=2,nb+1
|
374 scholl 1.6 if (pix(i).ne.0.0) then
375 bin(i)=bin(i)/pix(i)
376 c print*,"Z=",i,pix(i),bin(i)
377 endif
|
378 arta 1.1 enddo
379 bin(1)=2.0*bin(2)-bin(3)
380 c nb+2 bin from above loop only contains intensity values for r=rmax
381 bin(nb+2)=2.0*bin(nb+1)-bin(nb)
|
382 scholl 1.6 c print*,"X>",nb,bin(nb+2)
|
383 arta 1.1
384 c Make the Chebyschev fit
385 rminp=rmin-dr
386 rmaxp=rmax+dr
387 spflag=0
|
388 scholl 1.6
389 c print*,"3->",bin
390
|
391 scholl 1.4 c print*,"before: ",spflag,nb,rad,bin,cp,cpp
392 c print*,rminp," ",rmaxp," ",dr,c,nc
|
393 scholl 1.6 call chebft(rminp,rmaxp,c,nc,ifail)
394
|
395 scholl 1.4 c print*,"after: ",spflag,nb,rad,bin,cp,cpp
396 if(ifail.gt.0)then
397 c add by IS Oct5
398 ifail=24
399 return
400 endif
|
401 arta 1.1 call chder(rminp,rmaxp,c,cp,nc)
402 call chder(rminp,rmaxp,cp,cpp,nc)
403
404 c Do a least-squares calculation
405 lsq=0.0
406 ln=0.0
407 do 90 i=1,n
408 a=an(3,i)
409 if((a.le.alo).or.(a.gt.ahi)) goto 90
410 ln=ln+1.0
411 r=((an(1,i)-cmx)**2+(an(2,i)-cmy)**2)**0.5
412 lsq=lsq+(chebev(rminp,rmaxp,c,nc,r,ifail)-a)**2
413 90 enddo
|
414 scholl 1.4 if(ifail.gt.0) then
415 c add by IS Oct5
416 ifail=25
417 return
418 endif
|
419 arta 1.1 if(ln.eq.0.0)then
420 c This should not happen with normal data.
421 ifail=2
422 return
423 endif
424 dev=(lsq/ln)**0.5
425
426 c Cut points from data if they are "rem" sigma away from the fit
427 wbad=nbad ! save the number cut to compare to the next pass
428 nbad=0
429 cut=cut+1
430 c If "lmt" is exceeded, the parameter passed to output file "nbad" should be
431 c obviously large and "ifail"=3.
432 if(cut.gt.lmt)then
433 ifail=3
434 return
435 endif
436
437 vr=(dev*rem)**2
438 do 100 i=1,n
439 a=an(3,i)
440 arta 1.1 if((a.le.alo).or.(a.gt.ahi)) goto 100
441 r=((an(1,i)-cmx)**2+(an(2,i)-cmy)**2)**0.5
442 err=(a-chebev(rminp,rmaxp,c,nc,r,ifail))**2
443 if(err.ge.vr)then
444 nbad=nbad+1
445 an(3,i)=flag
446 endif
447 100 enddo
|
448 scholl 1.4 if(ifail.gt.0) then
449 c add by IS Oct5
450 ifail=26
451 return
452 endif
|
453 arta 1.1 ncut=ncut+nbad
454 c This criteria works ok for reasonable data.
455 if((wbad.eq.0).and.(nbad.eq.0)) goto 110 ! continue on
456 goto 60 ! take another cut
457
458 110 continue
459 c Compute the radial profile for "nang" angular bins
|
460 scholl 1.8 do j=1,jang
461 do i=1,jprf
462 prf(i,j)=0.0
463 enddo
464 enddo
|
465 arta 1.1 dprf=dr ! =(rmax-rmin)/float(nb)
466 do 118 i=1,n
467 a=an(3,i)
468 if((a.le.alo).or.(a.gt.ahi)) goto 118
469 x=an(1,i)
470 y=an(2,i)
471 r=((x-cmx)**2+(y-cmy)**2)**0.5
472 c this added to preshift limb distortion to get cleaner mean LDF
473 theta=atan2((y-cmy),(x-cmx))
474 if(theta.lt.0.0) theta=theta+tpi
475 ind=int(theta/dreg+1.0)
476 r = r - b0(ind)
477 c end of radius correction
478 jnd=int((r-rmin)/dprf+1.0)
479 if(jnd.lt.1)goto 118
480 if(jnd.gt.nprf)then
481 c Just in case the common radial range doesn't include all profiles
482 if(jnd.gt.nprf+2)then
483 c Need to consider making a new radial range...or there is problems in the data
484 ifail=4
485 return
486 arta 1.1 endif
487 jnd=nprf
488 endif
489 theta=atan2((y-cmy),(x-cmx))
490 if(theta.lt.0.0) theta=theta+tpi
491 ind=int(theta/dang+1.0)
492 if(ind.gt.nang)then
493 c This should not happen with normal data.
494 c ifail=5
495 ind=nang
496 c return
497 endif
498 lprf(jnd,ind)=lprf(jnd,ind)+a
499 prf(jnd,ind)=prf(jnd,ind)+1.0
500 118 enddo
501 do i=1,nang
502 do j=1,nprf
503 if(prf(j,i).gt.0.1) lprf(j,i)=lprf(j,i)/prf(j,i)
504 enddo
505 end do
506 c The last profile in the array is the average. "bin" & "rad" were offset
507 arta 1.1 c one array element to insure the Chebyschev fit do a good job near the edges
|
508 scholl 1.6 c print*,"jang=",jang
|
509 arta 1.1 do i=1,nprf
510 lprf(i,jang)=bin(i+1)
511 rprf(i)=rad(i+1)
|
512 scholl 1.6 c print*,"avg #:",i,lprf(i,jang),rprf(i)
|
513 arta 1.1 enddo
514 c Determine quantities per angular bin for alpha, beta calculation
515 do j=1,jreg
516 do i=1,6
517 expn(i,j)=0.0
518 enddo
519 enddo
520 do 120 i=1,n
521 a=an(3,i)
522 if((a.le.alo).or.(a.gt.ahi)) goto 120
523 x=an(1,i)
524 y=an(2,i)
525 r=((x-cmx)**2+(y-cmy)**2)**0.5
526 theta=atan2((y-cmy),(x-cmx))
527 if(theta.lt.0.0) theta=theta+tpi
528 ind=int(theta/dreg+1.0)
|
529 scholl 1.2 c IS Tue Sep 21 14:54:47 PDT 2010: switched the 2 following paragraphs
|
530 arta 1.1 if(ind.gt.nreg)then
531 ind=nang
532 c ifail=6
533 c return
534 endif
|
535 scholl 1.2 c correct radius bin
536 r = r - b0(ind)
537 c end correction
|
538 arta 1.1 D0=chebev(rminp,rmaxp,c,nc,R,ifail)
539 if(ifail.eq.11) then
|
540 scholl 1.4 c print*,"ifail=11!"
|
541 arta 1.1 ifail=0
542 goto 120
543 endif
544 D1=chebev(rminp,rmaxp,cp,nc,R,ifail)
545 D2=chebev(rminp,rmaxp,cpp,nc,R,ifail)
546 expn(1,ind)=expn(1,ind)+D0*a
547 expn(2,ind)=expn(2,ind)+D1*a
548 expn(3,ind)=expn(3,ind)+D2*a
549 expn(4,ind)=expn(4,ind)+D0*D1
550 expn(5,ind)=expn(5,ind)+D0*D2+D1**2
551 expn(6,ind)=expn(6,ind)+D0**2
552 120 enddo
|
553 scholl 1.4 if(ifail.gt.0)then
554 c add by IS Oct5
555 ifail=27
556 return
557 endif
|
558 arta 1.1 c Determine alpha (the LDF scale factor) and beta (the offset)
559 do i=1,nreg
560 thet(i)=float(i)*360.0/float(nreg)
561 DA=expn(1,i)
562 DB=expn(2,i)
563 DC=expn(3,i)
564 DD=expn(4,i)
565 DE=expn(5,i)
566 DF=expn(6,i)
567 DBB=DA*DE-DC*DF-DB*DD
568 if(DBB.eq.0.0)then
569 beta(i)=-10.0
570 alph(i)=-10.0
571 ifail=7
572 goto 130
573 c If too many points are cut, DBB (and DAA) can be zero. The flag this with
574 c ifail=7 but continue on. Alpha,betas with value -10 in a angular particular
575 c bin will show up in the output profiles. Any further software will have check
576 c for this unnatural phenomenon.
577 endif
578 beta(i)=(DA*DD-DB*DF)/DBB
579 arta 1.1 DBB=beta(i)
580 DAA=DD-DBB*DE
581 if(DAA.eq.0.0)then
582 alph(i)=-10.0
583 ifail=7
584 goto 130
585 endif
586 alph(i)=(DB-DBB*DC)/DAA-1.0
587 130 enddo
|
588 scholl 1.6 c print*,"end... ",cmx,cmy
|
589 arta 1.1 return
590 end
591
592
593 c Convert data into angle vs. radius, intensity, or I*r using trial center,
594 c cuts in radius, and intensity cuts
595 Subroutine Darr(anls,npts,ain,aout,cmx,cmy,rin,rout,thet,inte,n)
|
596 scholl 1.6 implicit none
597 integer jpt
598 parameter(jpt=8000000)
|
599 arta 1.1
600 integer i,n,npts
601 real anls(3,jpt),aout,ain,rout,rin,thet(1),inte(1)
602 real a,x,y,theta,r,tpi,cmx,cmy
603
604 tpi=6.283185307
605 n=0
606
607 do 10 i=1,npts
608 a=anls(3,i)
609 if((a.le.aout).or.(a.ge.ain)) goto 10
610 x=anls(1,i)
611 y=anls(2,i)
612 r=((x-cmx)**2+(y-cmy)**2)**0.5
613 if((r.lt.rin).or.(r.gt.rout)) goto 10
614 theta=atan2(cmx-x,cmy-y)
615 if(theta.lt.0.0) theta=theta+tpi
616 n=n+1
617 thet(n)=theta
618 c inte(n)=a*r
619 c inte(n)=a
620 arta 1.1 inte(n)=r
621 10 enddo
622
623 return
624 end
625
626
627 c Flags points to cut for a better fit
628 Subroutine Fitb(thet,inte,n,rem,chi,sig,cf,npts)
|
629 scholl 1.6 implicit none
|
630 arta 1.1 integer i,n,npts
631 real thet(1),inte(1),rem,chi,sig(1),cf(3),dev,err,tp
632
633 npts=0
634 dev=rem*(chi/float(n))**0.5
635 do i=1,n
636 tp=cf(1)*sin(thet(i))+cf(2)*cos(thet(i))+cf(3)
637 err=abs(tp-inte(i))
638 if(err.ge.dev)then
639 npts=npts+1
640 sig(i)=99999.0 ! sig(i) is normally 1.0
641 endif
642 enddo
643
644 return
645 end
646
647
648 c Computes the limb darkening function at any radius with a cubic spline
649 Real Function Ldf(r,ifail)
|
650 scholl 1.6 implicit none
651 integer jb
|
652 arta 1.1 parameter(jb=200)
653
654 integer nb,nb2,spflag,ifail
655 real r,bin(jb+2),rad(jb+2),snd(jb)
656
657 common /ldfcom/ spflag,nb,rad,bin
658
|
659 scholl 1.6 nb2=nb+2
|
660 arta 1.1 ldf=0.0
661 if(spflag.eq.0)then
662 call spline(rad,bin,nb2,2.0e30,2.0e30,snd)
663 spflag=1
664 endif
665 call splint(rad,bin,snd,nb2,r,ldf,ifail)
666
667 return
668 end
669
670
671 c Performs a Heap Sort (see Numerical Recipes) to find the mean of array "ra"
672 Subroutine Hpsort(n,ra)
|
673 scholl 1.6 implicit none
|
674 arta 1.1 integer i,ir,j,l,n
675 real ra(1),rra
676
677 if(n.lt.2) return
678
679 l=n/2+1
680 ir=n
681 10 continue
682
683 if(l.gt.1)then
684 l=l-1
685 rra=ra(l)
686 else
687 rra=ra(ir)
688 ra(ir)=ra(l)
689 ir=ir-1
690 if(ir.eq.1)then
691 ra(1)=rra
692 return
693 endif
694 endif
695 arta 1.1 i=l
696 j=l+1
697 20 if(j.le.ir)then
698 if(j.lt.ir)then
699 if(ra(j).lt.ra(j+1)) j=j+1
700 endif
701 if(rra.lt.ra(j))then
702 ra(i)=ra(j)
703 i=j
704 j=j+j
705 else
706 j=ir+1
707 endif
708 goto 20
709 endif
710 ra(i)=rra
711 goto 10
712
713 end
714
715
|
716 scholl 1.6 SUBROUTINE COVSRT(COVAR,NCVM,MA,LISTA,MFIT)
717 implicit none
718 integer ncvmp,ma,mfit
|
719 arta 1.1 parameter(ncvmp=3)
720 real COVAR(ncvmp,ncvmp)
|
721 scholl 1.6 real swap
722 integer LISTA(3),j,i,ncvm
723
|
724 arta 1.1 DO 12 J=1,MA-1
725 DO 11 I=J+1,MA
726 COVAR(I,J)=0.
727 11 CONTINUE
728 12 CONTINUE
729
730 DO 14 I=1,MFIT-1
731 DO 13 J=I+1,MFIT
732 IF(LISTA(J).GT.LISTA(I)) THEN
733 COVAR(LISTA(J),LISTA(I))=COVAR(I,J)
734 ELSE
735 COVAR(LISTA(I),LISTA(J))=COVAR(I,J)
736 ENDIF
737 13 CONTINUE
738 14 CONTINUE
739
740 SWAP=COVAR(1,1)
741 DO 15 J=1,MA
742 COVAR(1,J)=COVAR(J,J)
743 COVAR(J,J)=0.
744 15 CONTINUE
745 arta 1.1 COVAR(LISTA(1),LISTA(1))=SWAP
746 DO 16 J=2,MFIT
747 COVAR(LISTA(J),LISTA(J))=COVAR(1,J)
748 16 CONTINUE
749 DO 18 J=2,MA
750 DO 17 I=1,J-1
751 COVAR(I,J)=COVAR(J,I)
752 17 CONTINUE
753 18 CONTINUE
754
755 RETURN
756 END
757
758
|
759 scholl 1.6 SUBROUTINE GAUSSJ(A,N,NP,B,M,MP,ifail)
760 implicit none
761 integer nmax,npp,mpp,n,mp,ifail,m,np
|
762 arta 1.1 PARAMETER(NMAX=200,npp=3,mpp=1)
|
763 scholl 1.6 real B(npp,mpp),A(npp,npp),big,dum,pivinv
764 integer IPIV(NMAX),INDXR(NMAX),INDXC(NMAX),j,k,i,irow,icol,l,ll
|
765 arta 1.1
766 DO 11 J=1,N
767 IPIV(J)=0
768 11 CONTINUE
769
770 DO 22 I=1,N
771 BIG=0.
772 DO 13 J=1,N
773 IF(IPIV(J).NE.1)THEN
774 DO 12 K=1,N
775 IF (IPIV(K).EQ.0) THEN
776 IF (ABS(A(J,K)).GE.BIG) THEN
777 BIG=ABS(A(J,K))
778 IROW=J
779 ICOL=K
780 ENDIF
781 ELSE IF (IPIV(K).GT.1) THEN
782 ifail=9
783 return
784 ENDIF
785 12 CONTINUE
786 arta 1.1 ENDIF
787 13 CONTINUE
788 IPIV(ICOL)=IPIV(ICOL)+1
789
790 IF (IROW.NE.ICOL) THEN
791 DO 14 L=1,N
792 DUM=A(IROW,L)
793 A(IROW,L)=A(ICOL,L)
794 A(ICOL,L)=DUM
795 14 CONTINUE
796 DO 15 L=1,M
797 DUM=B(IROW,L)
798 B(IROW,L)=B(ICOL,L)
799 B(ICOL,L)=DUM
800 15 CONTINUE
801 ENDIF
802
803 INDXR(I)=IROW
804 INDXC(I)=ICOL
805 IF (A(ICOL,ICOL).EQ.0.) THEN
806 ifail=9
807 arta 1.1 return
808 endif
809 PIVINV=1./A(ICOL,ICOL)
810 A(ICOL,ICOL)=1.
811 DO 16 L=1,N
812 A(ICOL,L)=A(ICOL,L)*PIVINV
813 16 CONTINUE
814 DO 17 L=1,M
815 B(ICOL,L)=B(ICOL,L)*PIVINV
816 17 CONTINUE
817 DO 21 LL=1,N
818 IF(LL.NE.ICOL)THEN
819 DUM=A(LL,ICOL)
820 A(LL,ICOL)=0.
821 DO 18 L=1,N
822 A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
823 18 CONTINUE
824 DO 19 L=1,M
825 B(LL,L)=B(LL,L)-B(ICOL,L)*DUM
826 19 CONTINUE
827 ENDIF
828 arta 1.1 21 CONTINUE
829 22 CONTINUE
830 DO 24 L=N,1,-1
831 IF(INDXR(L).NE.INDXC(L))THEN
832 DO 23 K=1,N
833 DUM=A(K,INDXR(L))
834 A(K,INDXR(L))=A(K,INDXC(L))
835 A(K,INDXC(L))=DUM
836 23 CONTINUE
837 ENDIF
838 24 CONTINUE
839 RETURN
840 END
841
842
843 C FUNCTION FOR LFIT
844 SUBROUTINE FUNCS(X,AFUNC)
|
845 scholl 1.6 implicit none
846 REAL AFUNC(3) ,x
|
847 arta 1.1
848 AFUNC(1)=SIN(X)
849 AFUNC(2)=COS(X)
850 AFUNC(3)=1.0
851
852 RETURN
853 END
854
855
|
856 scholl 1.8 SUBROUTINE LFIT(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,COVAR,NCVM,CHISQ,
857 & FUNCS,ifail)
|
858 scholl 1.6 implicit none
859 integer mmax,ncvmp,jpt
860 PARAMETER (MMAX=3,ncvmp=3,jpt=8000000)
861 integer lista(3),kk,j,k,ihit,ifail,i,ncvm,ndata,mfit,ma
862 real X(jpt),Y(jpt),SIG(jpt),A(3),COVAR(ncvmp,ncvmp)
863 real BETA(MMAX),AFUNC(MMAX),CHISQ,SUM,WT,SIG2I,YM
864 external funcs
|
865 arta 1.1
866 KK=MFIT+1
867 DO 12 J=1,MA
868 IHIT=0
869 DO 11 K=1,MFIT
870 IF (LISTA(K).EQ.J) IHIT=IHIT+1
871 11 CONTINUE
872 IF (IHIT.EQ.0) THEN
873 LISTA(KK)=J
874 KK=KK+1
875 ELSE IF (IHIT.GT.1) THEN
876 ifail=8
877 return
878 ENDIF
879 12 CONTINUE
880 IF (KK.NE.(MA+1)) THEN
881 ifail=8
882 return
883 ENDIF
884
885 DO 14 J=1,MFIT
886 arta 1.1 DO 13 K=1,MFIT
887 COVAR(J,K)=0.
888 13 CONTINUE
889 BETA(J)=0.
890 14 CONTINUE
891
892 DO 18 I=1,NDATA
893 CALL FUNCS(X(I),AFUNC)
894 YM=Y(I)
895 IF(MFIT.LT.MA) THEN
896 DO 15 J=MFIT+1,MA
897 YM=YM-A(LISTA(J))*AFUNC(LISTA(J))
898 15 CONTINUE
899 ENDIF
900 SIG2I=1./SIG(I)**2
901 DO 17 J=1,MFIT
902 WT=AFUNC(LISTA(J))*SIG2I
903 DO 16 K=1,J
904 COVAR(J,K)=COVAR(J,K)+WT*AFUNC(LISTA(K))
905 16 CONTINUE
906 BETA(J)=BETA(J)+YM*WT
907 arta 1.1 17 CONTINUE
908 18 CONTINUE
909
910 IF (MFIT.GT.1) THEN
911 DO 21 J=2,MFIT
912 DO 19 K=1,J-1
913 COVAR(K,J)=COVAR(J,K)
914 19 CONTINUE
915 21 CONTINUE
916 ENDIF
917
918 CALL GAUSSJ(COVAR,MFIT,NCVM,BETA,1,1,ifail)
|
919 scholl 1.4 if(ifail.gt.0) then
920 c add by IS Oct5
921 ifail=30
922 return
923 endif
|
924 arta 1.1 DO 22 J=1,MFIT
925 A(LISTA(J))=BETA(J)
926 22 CONTINUE
927 CHISQ=0.
928 DO 24 I=1,NDATA
929 CALL FUNCS(X(I),AFUNC)
930 SUM=0.
931 DO 23 J=1,MA
932 SUM=SUM+A(J)*AFUNC(J)
933 23 CONTINUE
934 CHISQ=CHISQ+((Y(I)-SUM)/SIG(I))**2
935 24 CONTINUE
936
937 CALL COVSRT(COVAR,NCVM,MA,LISTA,MFIT)
938
|
939 scholl 1.6 RETURN
940 END
|
941 arta 1.1
942
|
943 scholl 1.6 SUBROUTINE SPLINE(X,Y,N,YP1,YPN,Y2)
944 implicit none
945 integer nmax,jb2
946 PARAMETER (NMAX=200,jb2=202)
947 real X(jb2),Y(jb2),Y2(jb2),U(NMAX),yp1,YPN,p,qn,un,k,sig
948 integer n,i
|
949 arta 1.1
950 IF (YP1.GT..99E30) THEN
951 Y2(1)=0.
952 U(1)=0.
953 ELSE
954 Y2(1)=-0.5
955 U(1)=(3./(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)
956 ENDIF
957
958 DO 11 I=2,N-1
959 SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1))
960 P=SIG*Y2(I-1)+2.
961 Y2(I)=(SIG-1.)/P
962 U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))
963 * /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P
964 11 CONTINUE
965
966 IF (YPN.GT..99E30) THEN
967 QN=0.
968 UN=0.
969 ELSE
970 arta 1.1 QN=0.5
971 UN=(3./(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))
972 ENDIF
973 Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.)
974 DO 12 K=N-1,1,-1
975 Y2(K)=Y2(K)*Y2(K+1)+U(K)
976 12 CONTINUE
977
|
978 scholl 1.6 RETURN
979 END
|
980 arta 1.1
981
|
982 scholl 1.6 SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y,ifail)
983 implicit none
984 integer jb2
985 PARAMETER (jb2=202)
986
987 real XA(jb2),YA(jb2),Y2A(jb2),x,y,a,b
988 integer k,klo,khi,h,ifail,n
|
989 arta 1.1
990 KLO=1
991 KHI=N
992 1 IF (KHI-KLO.GT.1) THEN
993 K=(KHI+KLO)/2
994 IF(XA(K).GT.X)THEN
995 KHI=K
996 ELSE
997 KLO=K
998 ENDIF
999 GOTO 1
1000 ENDIF
1001 H=XA(KHI)-XA(KLO)
|
1002 scholl 1.6 IF(H.EQ.0.)then
1003 ifail=10
1004 c print*,"ifail=10"
1005 return
|
1006 arta 1.1 endif
1007 A=(XA(KHI)-X)/H
1008 B=(X-XA(KLO))/H
1009 Y=A*YA(KLO)+B*YA(KHI)+
1010 * ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.
1011
|
1012 scholl 1.6 RETURN
1013 END
|
1014 arta 1.1
1015
1016 c Chebyschev Approximation to function expressed in func on interval
1017 c [a,b]. Numerical Recipes pp 184-190.
1018
1019
1020 c Compute coefficients c(1:n)
|
1021 scholl 1.6 Subroutine CHEBFT(a,b,c,n,ifail)
1022 implicit none
1023 integer nmax
1024 real pi
|
1025 arta 1.1 parameter(nmax=200,pi=3.141592653589793D0)
1026 integer n,j,k,ifail
|
1027 scholl 1.6 real a,b,c(nmax),bma,bpa,fac,y,f(nmax),ldf
|
1028 arta 1.1 double precision sum
|
1029 scholl 1.6 c external func
|
1030 arta 1.1
1031 bma=0.5*(b-a)
1032 bpa=0.5*(b+a)
1033 do k=1,n
1034 y=cos(pi*(k-0.5)/n)
|
1035 scholl 1.6 f(k)=ldf(y*bma+bpa,ifail)
|
1036 arta 1.1 enddo
|
1037 scholl 1.4 if(ifail.gt.0) then
1038 c add by IS Oct5
1039 ifail=31
|
1040 scholl 1.6 c print*,"ifail=31"
|
1041 scholl 1.4 return
1042 endif
|
1043 arta 1.1 fac=2.0/n
1044 do j=1,n
1045 sum=0.D0
1046 do k=1,n
1047 sum=sum+f(k)*cos((pi*(j-1))*((k-0.5D0)/n))
1048 enddo
1049 c(j)=fac*sum
1050 enddo
1051
1052 return
1053 end
1054
1055
1056 c Evaluate the approximation to f(x)
1057 Function CHEBEV(a,b,c,m,x,ifail)
|
1058 scholl 1.6 implicit none
|
1059 arta 1.1
|
1060 scholl 1.6 integer m,j,ifail,jb
1061 parameter(jb=200)
1062 real chebev,a,b,x,c(jb),d,dd,sv,y,y2
|
1063 arta 1.1
1064 if((x-a)*(x-b).gt.0.0) then
1065 c CHEBEV out of range
1066 ifail=11
1067 return
1068 endif
1069 d=0.0
1070 dd=0.0
1071 y=(2.0*x-a-b)/(b-a)
1072 y2=2.0*y
1073 do j=m,2,-1
1074 sv=d
1075 d=y2*d-dd+c(j)
1076 dd=sv
1077 enddo
1078 chebev=y*d-dd+0.5*c(1)
1079
1080 return
1081 end
1082
1083
1084 arta 1.1 c Compute the coefficients of the derivative of the function whose
1085 c coefficients are c(n)
1086 Subroutine CHDER(a,b,c,cder,n)
|
1087 scholl 1.6 implicit none
1088 integer jc
1089 parameter(jc=200)
|
1090 arta 1.1 integer n,j
|
1091 scholl 1.6 real a,b,c(jc),cder(jc),con
|
1092 arta 1.1
1093 cder(n)=0.0
1094 cder(n-1)=2.0*(n-1)*c(n)
1095 if(n.ge.3)then
1096 do j=n-2,1,-1
1097 cder(j)=cder(j+2)+2.0*j*c(j+1)
1098 enddo
1099 endif
1100 con=2.0/(b-a)
1101 do j=1,n
1102 cder(j)=cder(j)*con
1103 enddo
1104
1105 return
1106 end
1107
|
Return to
limb.f
CVS log
Up to [Development] / JSOC / proj / limbfit / apps