GeneralBrokenLines  Rev:70
 All Classes Files Functions Variables Pages
gbltraj.f90
Go to the documentation of this file.
1 !
2 
3 ! Code converted using TO_F90 by Alan Miller
4 ! Date: 2012-03-03 Time: 16:59:41
5 
6 ! C. Kleinwort, DESY-FH1, www.terascale.de
7 ! March 2011, Analysis Centre: Statistics Tools Group
8 
82 
84 MODULE gbltraj
85  IMPLICIT NONE
86  SAVE
87  INTEGER, PARAMETER :: maxFitPar=250 !< max. number of fit parameters
88  INTEGER, PARAMETER :: maxFitSymMatSize=(maxFitPar*maxFitPar+maxFitPar)/2
89  INTEGER, PARAMETER :: maxBandWidth=10 !< max. band width
90  INTEGER, PARAMETER :: maxBorderSize=10 !< max. border size
91  INTEGER, PARAMETER :: maxTrackPar=maxBorderSize+5 !< max. number of local track parameters
92  INTEGER, PARAMETER :: maxTrackSymMatSize=(maxTrackPar*maxTrackPar+maxTrackPar)/2
93  INTEGER, PARAMETER :: maxPoints=1000 !< max. points
94  INTEGER, PARAMETER :: sizeDataBuffer=10*maxPoints !< max. integrated size of data blocks
95  !
96  INTEGER :: trajLevel !< level of (preparation of) trajectory (fit)
97  INTEGER :: printLevel !< print level
98  INTEGER :: offsetDimension !< dimension of offsets
99  !! (2: 2D offsets for track in 3D, 1: 1D offsets for track in 2D)
100  INTEGER :: numPoints !< number of points
101  INTEGER :: numMeas !< number (of points) with measurement
102  INTEGER :: numScat !< number (of points) with scatterer
103  INTEGER :: numOffsets !< number (of points) with offsets (as fit parameter)
104  INTEGER :: numKinks !< number of (multiple scattering) kinks
105  INTEGER :: lastOffset !< last point with offset
106  INTEGER :: lastScatterer !< last point with scatterer
107  INTEGER :: numFitPar !< number of fit parameters
108  INTEGER :: numCurv !< flag (0/1) for usage of Q/P as fit parameter
109  INTEGER :: bandSize !< band size
110  INTEGER :: numAddLocPar !< number of additional local track parameters
111  INTEGER :: extSeedPoint !< point with external seed (or 0)
112  INTEGER :: firstActiveCoord !< first active offset coordinate (0: u_1, 1: u_2)
113  INTEGER :: lastDownWeightMethod !< last used down-weighting method
114  DOUBLE PRECISION, DIMENSION(maxFitPar) :: vecb !< right hand side of linear equation system (A*x=b)
115  DOUBLE PRECISION, DIMENSION(maxFitSymMatSize) :: matA !< (sym.) matrix of linear equation system (A*x=b)
116  DOUBLE PRECISION, DIMENSION(maxTrackSymMatSize) :: extSeedMat !< external seed (precision matrix)
117  DOUBLE PRECISION, DIMENSION(maxTrackPar*maxTrackPar) :: jacTrackToFit !< jacobian for transformation fit to track parameter
118 
119  TYPE gblpoint
120  INTEGER :: offset !< offset at point (or from previous point)
121  INTEGER :: indMeas !< measurement flag
122  INTEGER :: indScat !< scatterer flag
123  INTEGER :: indKink !< kink index
124  DOUBLE PRECISION, DIMENSION(5,5) :: jacPointToPoint !< point-to-point jacobian (from previous point)
125  DOUBLE PRECISION, DIMENSION(2,2) :: prevJ !< matrix 'J', offset part of jacobian to previous point with offset
126  DOUBLE PRECISION, DIMENSION(2,2) :: prevS !< matrix 'S', slope part of jacobian to previous point with offset
127  DOUBLE PRECISION, DIMENSION(2) :: prevd !< vector 'd', Q/P part of jacobian to previous point with offset
128  DOUBLE PRECISION, DIMENSION(2,2) :: nextJ !< matrix 'J', offset part of jacobian to next point with offset
129  DOUBLE PRECISION, DIMENSION(2,2) :: nextS !< matrix 'S', slope part of jacobian to next point with offset
130  DOUBLE PRECISION, DIMENSION(2) :: nextd !< vector 'd', Q/P part of jacobian to next point with offset
131  DOUBLE PRECISION, DIMENSION(2,2) :: matProj !< projection matrix from measurement to local coordinate system
132  INTEGER, DIMENSION(2) :: indLocal !< index of local derivatives
133  INTEGER, DIMENSION(2) :: indGlobal !< index of global derivatives
134  REAL, DIMENSION(2) :: measValue !< residual (for measurement)
135  REAL, DIMENSION(2) :: measPrec !< precision (for measurement, 1/sigma**2)
136  REAL, DIMENSION(2) :: scatValue !< residual (for kink)
137  REAL, DIMENSION(2) :: scatPrec !< precision (for kink, 1/sigma**2)
138  END TYPE gblpoint
139  TYPE(gblpoint), DIMENSION(maxPoints) :: points !< list of GblPoints
140 
141  INTEGER :: numDataBlocks !< number of data blocks with measurements or kinks
142  INTEGER :: integDataSize !< integrated size of data blocks
143  INTEGER :: maxDataSize !< max. integrated size
144  INTEGER, DIMENSION(sizeDataBuffer) :: intData !< integer part of data blocks (lower part, 1..mdat)
145  !! or local/global derivatives (upper part, mdat+1..mxdat)
146  REAL, DIMENSION(sizeDataBuffer) :: floatData !< float part of data blocks or local/global derivatives
147 
148  INTEGER :: maxUsedFitPar = 0 !< max. number of fit parameters used
149  INTEGER :: lastNumFitPar = 0 !< number of fit parameters from last fit
150  INTEGER :: lastBorderSize = 0 !< border size from last fit
151  INTEGER :: lastBandWidth = 0 !< band width from last fit
152 !
153 CONTAINS
154 
166  SUBROUTINE gblfit(cdw,mrank,np,ndf,chi2,wls)
167  IMPLICIT NONE
168  REAL :: als
169  INTEGER :: idwm
170  INTEGER :: inv
171  INTEGER :: iter
172  INTEGER :: lcdw
173  INTEGER :: nrank
174  REAL :: swgt
175 
176  CHARACTER (LEN=*), INTENT(IN) :: cdw
177  INTEGER, INTENT(OUT) :: mrank
178  INTEGER, INTENT(OUT) :: np
179  INTEGER, INTENT(OUT) :: ndf
180  REAL, INTENT(OUT) :: chi2
181  REAL, INTENT(OUT) :: wls
182 
183  mrank=0
184  np=0
185  ndf=-1
186  chi2=0.
187  wls=0.
188  lcdw=len_trim(cdw)
189  iter=0
190 
191  CALL gblprp
192  IF (trajlevel < 3) return
193 
194  iterate: DO
195  CALL gblmat
196  mrank=numcurv+2*offsetdimension ! rank of trajectory
197  np =5+numaddlocpar ! number of track parameters
198 
199  IF (printlevel > 1) print *, ' GBLFIT: NPAR, NBND, NBDR ', numfitpar, bandsize, numcurv+numaddlocpar
200  inv=1 ! calc band part of covariance matrix
201  CALL sqmibb(mata,vecb,numfitpar,numcurv+numaddlocpar,bandsize,inv,nrank)
202  IF (nrank /= numfitpar) THEN ! rank deficit
203  IF (printlevel > 1) print *, ' GBLFIT: rank deficit ', nrank, numfitpar
204  mrank=-mrank
205  chi2=0.
206  ndf=0
207  als=0.
208  return
209  END IF
210 
211  iter=iter+1
212  idwm=0 ! down weighting method (Tukey, Huber, Cauchy)
213  IF (iter <= lcdw) idwm=(index('TtHhCc',cdw(iter:iter))+1)/2
214 
215  CALL gblch2(idwm,chi2,swgt)
216  ndf=numdatablocks-numfitpar
217  ! iterate ?
218  IF (idwm <= 0) exit
219  END DO iterate
220 
221  wls=float(numdatablocks)-swgt
222  trajlevel=4 ! fitted
223 
224  return
225  END SUBROUTINE gblfit
226 
229 
230  SUBROUTINE gblprp
231  IMPLICIT NONE
232  INTEGER :: i
233  INTEGER :: ijsym
234  INTEGER :: j
235 
236  ! index in symmetric matrix
237  ijsym(i,j)=min(i,j)+(max(i,j)*max(i,j)-max(i,j))/2
238 
239  IF (trajlevel == 1) CALL gblcjc ! calc jacobians to prev/next offset
240  IF (trajlevel /= 2) return ! no points with jacobians added
241  IF (numoffsets < 2) return ! too few parameters
242  ! number of parameters
243  numaddlocpar=min(numaddlocpar,maxbordersize-numcurv) ! limit number of (external) local par.
244  numfitpar=numoffsets*offsetdimension+numcurv+numaddlocpar ! number of parameters
245  IF (numfitpar > maxfitpar) return ! too many parameters
246 
247  CALL gbldat
248 
249  trajlevel=3 ! prepared
250 
251  return
252  END SUBROUTINE gblprp
253 
256  SUBROUTINE gblcjc
257  IMPLICIT NONE
258 
259  INTEGER :: ipoint
260  INTEGER :: nstep
261  INTEGER :: i
262  INTEGER :: j
263  INTEGER :: lastpoint
264 
265  DOUBLE PRECISION :: ajaci(5,5)
266  DOUBLE PRECISION :: ajacn(5,5)
267  DOUBLE PRECISION :: ajact(5,5)
268 
269  ! forward propagation (all)
270 
271  nstep=0
272  lastpoint=1 ! last point with offset
273 
274  DO ipoint=2,numpoints
275  ! full jacobian (to next point)
276  ajacn=points(ipoint)%jacPointToPoint
277 
278  IF (nstep == 0) THEN ! start integrated jacobian
279  ajaci=ajacn
280  ELSE ! update integrated jacobian
281  forall(i=1:5,j=1:5) ajact(i,j)=sum(ajacn(i,:)*ajaci(:,j))
282  ajaci=ajact ! MATMUL(ajacn,ajaci) is slower
283  END IF
284  nstep=nstep+1
285  CALL gblinu(ajaci,ajact)
286  ! IPNT -> PREV
287  points(ipoint)%prevJ=ajact(4:5,4:5) ! J-
288  points(ipoint)%prevS=ajact(4:5,2:3) ! S-
289  points(ipoint)%prevd=ajact(4:5,1) ! d-
290 
291  IF (points(ipoint)%offset > 0) THEN
292  ! LPNT -> NEXT
293  points(lastpoint)%nextJ=ajaci(4:5,4:5) ! J+
294  points(lastpoint)%nextS=ajaci(4:5,2:3) ! S+
295  points(lastpoint)%nextd=ajaci(4:5,1) ! d+
296 
297  nstep=0 ! restart integration
298  lastpoint=ipoint
299  END IF
300  END DO
301 
302  ! backward propagation (without scatterers)
303 
304  nstep=0
305 
306  DO ipoint=numpoints-1,1,-1
307 
308  IF (points(ipoint)%offset > 0) THEN ! skip offsets
309  nstep=0
310  cycle
311  END IF
312  ! full jacobian (to next point)
313  ajacn=points(ipoint+1)%jacPointToPoint
314 
315  IF (nstep == 0) THEN ! start integrated jacobian
316  ajaci=ajacn
317  ELSE ! update integrated jacobian
318  forall(i=1:5,j=1:5) ajact(i,j)=sum(ajaci(i,:)*ajacn(:,j))
319  ajaci=ajact ! MATMUL(ajaci,ajacn) is slower
320  END IF
321  nstep=nstep+1
322  ! IPNT -> NEXT
323  points(ipoint)%nextJ=ajaci(4:5,4:5) ! J+
324  points(ipoint)%nextS=ajaci(4:5,2:3) ! S+
325  points(ipoint)%nextd=ajaci(4:5,1) ! d+
326  END DO
327 
328  trajlevel=2
329  return
330  END SUBROUTINE gblcjc
331 
337 
338  SUBROUTINE gblini(lprnt,icoord)
339  IMPLICIT NONE
340  INTEGER :: ifirst
341 
342  INTEGER, INTENT(IN) :: lprnt
343  INTEGER, OPTIONAL, INTENT(IN) :: icoord
344 
345  DATA ifirst / 1 /
346 
347  offsetdimension=2 ! 2D offsets for track in 3D
348  firstactivecoord=0
349 
350  IF (present(icoord)) THEN
351  IF (icoord > 0) THEN
352  offsetdimension=1 ! 1D offsets for track in 2D
353  firstactivecoord=0
354  IF (icoord > 1) firstactivecoord=1
355  END IF
356  END IF
357 
358  trajlevel=0
359  printlevel=lprnt
360  numcurv=1 ! with Q/P
361  numaddlocpar=0
362  numpoints=0
363  lastoffset=0
364  nummeas=0
365  numscat=0
366  lastscatterer=0
367  numoffsets=0
368  numkinks=0
369  maxdatasize=sizedatabuffer
370  extseedpoint=0 ! external seed parameter offset
371 
372  IF (ifirst == 1) THEN
373  ifirst=0
374  IF (printlevel > 0) print *, ' GBL $Rev: 67 $'
375  END IF
376 
377  return
378  END SUBROUTINE gblini
379 
385 
386  SUBROUTINE gbladp(ajac,ipoint)
387  IMPLICIT NONE
388 
389  DOUBLE PRECISION, INTENT(IN) :: ajac(5,5)
390  INTEGER, INTENT(OUT) :: ipoint
391 
392  ipoint=0
393  IF (numpoints >= maxpoints) THEN
394  IF (printlevel > 0) print *, ' GBLADP: too many points, ignoring new point '
395  return
396  END IF
397  numpoints=numpoints+1
398  points(numpoints)%offset=0
399  points(numpoints)%indMeas=0
400  points(numpoints)%indScat=0
401  points(numpoints)%indKink=0
402  points(numpoints)%indLocal=0
403  points(numpoints)%indGlobal=0
404  ipoint=numpoints
405 
406  points(numpoints)%jacPointToPoint=ajac
407 
408  points(numpoints)%prevJ=0.0d0
409  points(numpoints)%prevS=0.0d0
410  points(numpoints)%prevd=0.0d0
411  points(numpoints)%nextJ=0.0d0
412  points(numpoints)%nextS=0.0d0
413  points(numpoints)%nextd=0.0d0
414 
415  IF (lastoffset > 0) THEN
416  IF (points(lastoffset)%indScat <= 0.AND.numoffsets > 1) THEN
417  numoffsets=numoffsets-1
418  points(lastoffset)%offset=-points(lastoffset)%offset
419  END IF
420  END IF
421  numoffsets=numoffsets+1
422  points(numpoints)%offset=numoffsets
423  lastoffset=numpoints
424  IF (numoffsets > 2.AND.lastscatterer > 0) THEN
425  numkinks=numkinks+1
426  points(lastscatterer)%indKink=numkinks
427  lastscatterer=0
428  END IF
429 
430  trajlevel=1 ! adding points
431  return
432  END SUBROUTINE gbladp
433 
441 
442  SUBROUTINE gbladm(proj,res,prec)
443  IMPLICIT NONE
444 
445  DOUBLE PRECISION, INTENT(IN) :: proj(2,2)
446  REAL, INTENT(IN) :: res(2)
447  REAL, INTENT(IN) :: prec(2)
448 
449  IF (numpoints <= 0) THEN
450  IF (printlevel > 0) print *, ' GBLADM: no point defined'
451  return
452  END IF
453  IF (nummeas+numscat >= maxpoints) THEN
454  IF (printlevel > 0) print *, &
455  ' GBLADM: too many measurement+scatterers ', nummeas+numscat
456  return
457  END IF
458  IF (points(numpoints)%indMeas <= 0) THEN
459  nummeas=nummeas+1
460  points(numpoints)%indMeas=1
461  points(numpoints)%matProj=proj
462  points(numpoints)%measValue=res
463  points(numpoints)%measPrec=prec
464  ELSE
465  IF (printlevel > 0) print *, &
466  ' GBLADM: measurement already defined for point ', numpoints
467  END IF
468 
469  return
470  END SUBROUTINE gbladm
471 
477 
478  SUBROUTINE gbladl(nder,der,iret)
479  IMPLICIT NONE
480  INTEGER :: im
481  INTEGER :: k
482  INTEGER :: mdat0
483 
484  INTEGER, INTENT(IN) :: nder
485  REAL, INTENT(IN) :: der(2,nder)
486  INTEGER, INTENT(OUT) :: iret
487 
488  iret=0
489  IF (numpoints <= 0) THEN
490  IF (printlevel > 0) print *, ' GBLADL: no point defined'
491  return
492  END IF
493 
494  DO im=1,2
495 
496  IF (points(numpoints)%indLocal(im) == 0) THEN ! local derivs not yet defined
497 
498  IF (integdatasize > maxdatasize-nder) return ! enough space left
499  mdat0=maxdatasize
500  DO k=nder,1,-1
501  IF (der(im,k) /= 0.0) THEN
502  numaddlocpar=max(numaddlocpar,k)
503  intdata(maxdatasize)=k
504  floatdata(maxdatasize)=der(im,k)
505  maxdatasize=maxdatasize-1
506  END IF
507  END DO
508 
509  points(numpoints)%indLocal(im)=maxdatasize
510  intdata(maxdatasize)=mdat0-maxdatasize
511  floatdata(maxdatasize)=0.0
512  maxdatasize=maxdatasize-1
513  iret=iret+mdat0-maxdatasize
514 
515  ELSE
516  IF (printlevel > 0) print *, &
517  ' GBLADL: local derivatives already defined for point ', numpoints
518  END IF
519 
520  END DO
521 
522  return
523  END SUBROUTINE gbladl
524 
531 
532  SUBROUTINE gbladg(nder,lder,der,iret)
533  IMPLICIT NONE
534  INTEGER :: im
535  INTEGER :: k
536  INTEGER :: mdat0
537 
538  INTEGER, INTENT(IN) :: nder
539  INTEGER, INTENT(IN) :: lder(nder)
540  REAL, INTENT(IN) :: der(2,nder)
541  INTEGER, INTENT(OUT) :: iret
542 
543  iret=0
544  IF (numpoints <= 0) THEN
545  IF (printlevel > 0) print *, ' GBLADG: no point defined'
546  return
547  END IF
548 
549  DO im=1,2
550 
551  IF (points(numpoints)%indGlobal(im) == 0) THEN ! local derivs not yet defined
552 
553  IF (integdatasize > maxdatasize-nder) return ! enough space left
554  mdat0=maxdatasize
555  DO k=nder,1,-1
556  IF (der(im,k) /= 0.0) THEN
557  intdata(maxdatasize)=lder(k)
558  floatdata(maxdatasize)=der(im,k)
559  maxdatasize=maxdatasize-1
560  END IF
561  END DO
562 
563  points(numpoints)%indGlobal(im)=maxdatasize
564  intdata(maxdatasize)=mdat0-maxdatasize
565  floatdata(maxdatasize)=0.0
566  maxdatasize=maxdatasize-1
567  iret=iret+mdat0-maxdatasize
568 
569  ELSE
570  IF (printlevel > 0) print *, &
571  ' GBLADG: global derivatives already defined for point ', numpoints
572  END IF
573 
574  END DO
575 
576  return
577  END SUBROUTINE gbladg
578 
586 
587  SUBROUTINE gblads(res,prec)
588  IMPLICIT NONE
589  INTEGER :: jscat
590 
591  REAL, INTENT(IN) :: res(2)
592  REAL, INTENT(IN) :: prec(2)
593 
594  IF (numpoints <= 0) THEN
595  IF (printlevel > 0) print *, ' GBLADS: no point defined'
596  return
597  END IF
598  IF (nummeas+numscat >= maxpoints) THEN
599  IF (printlevel > 0) print *, &
600  ' GBLADS: too many measurement+scatterers ', nummeas+numscat
601  return
602  END IF
603  IF (prec(1) <= 0.0.OR.prec(2) <= 0.0) THEN
604  IF (printlevel > 0) print *, ' GBLADS: invalid scattering precision ', prec
605  return
606  END IF
607 
608  IF (points(numpoints)%indScat <= 0) THEN
609  jscat=maxpoints-numscat
610  numscat=numscat+1
611  points(numpoints)%indScat=1
612  lastscatterer=numpoints
613  points(numpoints)%scatValue=res
614  points(numpoints)%scatPrec=prec
615  ELSE
616  IF (printlevel > 0) print *, ' GBLADM: scatterer already defined for point ', numpoints
617  END IF
618 
619  return
620  END SUBROUTINE gblads
621 
623 
624  SUBROUTINE gbldmp
625  IMPLICIT NONE
626  INTEGER :: i
627 
628  IF (numpoints <= 0) THEN
629  print *, ' GBLDMP no trajectory defined '
630  return
631  END IF
632 
633  print *
634  print *, ' GBLDMP trajectory definition '
635  print *, '-------------------------------------'
636  print *, ' number of points ', numpoints
637  print *, ' number of offsets ', numoffsets
638  print *, ' number of scatterers ', numscat
639  print *, ' number of measurements ', nummeas
640  print *, ' number of kinks ', numkinks
641  IF (numaddlocpar > 0) print *, ' number of local par. ', numaddlocpar
642  print *
643 
644  DO i=1,numpoints
645  print *, ' Point ', i
646  IF (points(i)%offset > 0) print *, ' Offset ', points(i)%offset
647  IF (points(i)%indScat > 0) print *, ' Scatterer ', points(i)%indScat
648  IF (points(i)%indMeas > 0) THEN
649  IF (points(i)%offset > 0) THEN
650  print *, ' Measurement using offset ',points(i)%offset
651  ELSE
652  print *, ' Measurement using offsets ', &
653  -points(i)%offset-1,' ,',-points(i)%offset
654  END IF
655  END IF
656  IF (points(i)%indKink > 0) print *, ' Kink ', points(i)%indKink, &
657  ' using offsets ',points(i)%offset-1,'..',points(i)%offset+1
658  IF (points(i)%indMeas <= 0.AND.points(i)%indScat <= 0) THEN
659  IF (points(i)%offset > 0) THEN
660  print *, ' Prediction ', ' using offset ',points(i)%offset
661 
662  ELSE
663  print *, ' Prediction ', &
664  ' using offsets ',-points(i)%offset-1,' ,',-points(i)%offset
665  END IF
666  END IF
667  END DO
668 
669  return
670  END SUBROUTINE gbldmp
671 
673 
674  SUBROUTINE gbldat
675  IMPLICIT NONE
676  INTEGER :: i
677  INTEGER :: idx
678  INTEGER :: ilcl
679  INTEGER :: ioff1
680  INTEGER :: ip0
681  INTEGER :: ipar0
682  INTEGER :: ipmn
683  INTEGER :: ipoff
684  INTEGER :: j
685  INTEGER :: jmeas
686  INTEGER :: joff
687  INTEGER :: jpnt
688  INTEGER :: jscat
689  INTEGER :: k
690  INTEGER :: kdat
691  INTEGER :: kdat0
692  INTEGER :: l
693  INTEGER :: l0
694  INTEGER :: ldat0
695  INTEGER :: ltem
696  INTEGER :: m
697  INTEGER :: nb
698  INTEGER :: nlcl
699  INTEGER :: nplc
700  INTEGER :: npsd
701 
702  DOUBLE PRECISION :: tc2l(2,2)
703  DOUBLE PRECISION :: wp(2,2)
704  DOUBLE PRECISION :: wn(2,2)
705  DOUBLE PRECISION :: wjp(2,2)
706  DOUBLE PRECISION :: wjn(2,2)
707  DOUBLE PRECISION :: dwp(2)
708  DOUBLE PRECISION :: dwn(2)
709  DOUBLE PRECISION :: wjs(2,2)
710  DOUBLE PRECISION :: wji(2,2)
711  DOUBLE PRECISION :: pn(2,2)
712  DOUBLE PRECISION :: wpp(2,2)
713  DOUBLE PRECISION :: wpn(2,2)
714  DOUBLE PRECISION :: dws(2)
715  DOUBLE PRECISION :: dps(2)
716  DOUBLE PRECISION :: det
717  DOUBLE PRECISION :: diag(maxtrackpar)
718  DOUBLE PRECISION :: eigen(maxtrackpar*maxtrackpar)
719  DOUBLE PRECISION :: work(maxtrackpar)
720  DOUBLE PRECISION :: seed(maxtrackpar*maxtrackpar)
721  INTEGER :: iwork(maxtrackpar)
722 
723  l0=firstactivecoord
724 
725  lastdownweightmethod=0 ! start without downweighting
726  ! loop over points
727  numdatablocks=0 ! number of data items
728  integdatasize=0 ! length of data
729  ipmn=numfitpar+1 ! minimum parameter with non zero derivatives
730  DO jpnt=1,numpoints
731 
732  joff=points(jpnt)%offset ! offset at point
733  ipar0=offsetdimension*(joff-1)+numcurv+numaddlocpar ! offset in parameter number
734  ! kink at point ?
735  IF (points(jpnt)%indKink > 0) THEN
736  jscat=points(jpnt)%indScat
737 
738  CALL gblder(jpnt,-1,wp,wjp,dwp)
739  CALL gblder(jpnt, 2,wn,wjn,dwn)
740  dws=dwp+dwn ! W- * d- + W+ * d+
741  wjs=wjp+wjn ! W- * J- + W+ * J+
742 
743  DO m=1,offsetdimension ! pseudo measurements
744  kdat=maxdatasize-(numcurv+3*offsetdimension)
745  IF (integdatasize+3 > kdat) return
746  kdat0=kdat
747 
748  ldat0=integdatasize
749  numdatablocks=numdatablocks +1
750  intdata(integdatasize+1)=3 ! length of item
751  intdata(integdatasize+2)=jpnt ! point
752  intdata(integdatasize+3)=-m ! pseudo measurement
753  floatdata(integdatasize+1)=points(jpnt)%scatValue(m) ! (start) value
754  floatdata(integdatasize+2)=points(jpnt)%scatPrec(m) ! precision (1/error**2)
755  floatdata(integdatasize+3)=1.0 ! (additional) weight
756  integdatasize=integdatasize+3
757 
758  intdata(kdat+1)= 1
759  floatdata(kdat+1)= -sngl(dws(m))
760  kdat=kdat+numcurv
761  DO l=1,offsetdimension
762  idx=kdat+l
763  intdata(idx)= ipar0+l-offsetdimension
764  floatdata(idx)= sngl( wp(m,l+l0))
765  idx=idx+offsetdimension
766  intdata(idx)= ipar0+l
767  floatdata(idx)= sngl(-wjs(m,l+l0))
768  idx=idx+offsetdimension
769  intdata(idx)= ipar0+l+offsetdimension
770  floatdata(idx)= sngl( wn(m,l+l0))
771  END DO
772  kdat=kdat+3*offsetdimension
773 
774  DO k=kdat0+1,kdat
775  IF (floatdata(k) /= 0.0) THEN ! copy non zero derivatives
776  integdatasize=integdatasize+1
777  intdata(integdatasize)=intdata(k)
778  floatdata(integdatasize)=floatdata(k)
779  ipmn=min(ipmn,intdata(integdatasize))
780  END IF
781  END DO
782  intdata(ldat0+1)=integdatasize-ldat0
783  END DO
784 
785  END IF
786  ! measurement at point ?
787  jmeas=points(jpnt)%indMeas ! measurement at point
788  IF (jmeas > 0) THEN
789  tc2l=points(jpnt)%matProj ! P
790  ipar0=offsetdimension*(joff-1)+numcurv+numaddlocpar ! offset in parameter number
791  IF (joff < 0) THEN ! need interpolation
792  ipar0=offsetdimension*(-joff-2)+numcurv+numaddlocpar ! offset in parameter number
793  CALL gblder(jpnt,-1,wp,wjp,dwp)
794  CALL gblder(jpnt, 2,wn,wjn,dwn)
795  dws=dwp+dwn ! W- * d- + W+ * d+
796  wjs=wjp+wjn ! W- * J- + W+ * J+
797  det=wjs(1,1)*wjs(2,2)-wjs(1,2)*wjs(2,1) ! (W- * J- + W+ * J+)^-1 (=N)
798  wji(1,1)= wjs(2,2)/det;wji(1,2)=-wjs(1,2)/det
799  wji(2,1)=-wjs(2,1)/det;wji(2,2)= wjs(1,1)/det
800  forall(i=1:2,j=1:2) pn(i,j)=sum(tc2l(i,:)*wji(:,j)) ! P * N
801  forall(i=1:2,j=1:2) wpp(i,j)=sum(pn(i,:)*wp(:,j)) ! P * N * W-
802  forall(i=1:2,j=1:2) wpn(i,j)=sum(pn(i,:)*wn(:,j)) ! P * N * W+
803  forall(i=1:2) dps(i)=sum(dws(:)*pn(i,:)) ! P * N * (W+ * d+ + W- * d-)
804 
805  END IF
806 
807  DO m=1,2
808  IF (points(jpnt)%measPrec(m) <= 0.0) cycle ! no precision ?
809  kdat=maxdatasize-(numcurv+2*offsetdimension)
810  IF (integdatasize+3 > kdat) return
811  kdat0=kdat
812 
813  ldat0=integdatasize
814  numdatablocks=numdatablocks+1
815  intdata(integdatasize+1)=3 ! length of item
816  intdata(integdatasize+2)=jpnt ! point
817  intdata(integdatasize+3)=m ! measurement
818  floatdata(integdatasize+1)=points(jpnt)%measValue(m) ! value
819  floatdata(integdatasize+2)=points(jpnt)%measPrec(m) ! precision (1/error**2)
820  floatdata(integdatasize+3)=1.0 ! (additional) weight
821  integdatasize=integdatasize+3
822 
823  IF (joff > 0) THEN ! measurement at offset
824  DO l=1,offsetdimension
825  intdata(kdat+1)= ipar0+l
826  floatdata(kdat+1)= sngl(tc2l(m,l+l0))
827  kdat=kdat+1
828  END DO
829  ELSE ! interpolation between offsets
830 
831  intdata(kdat+1)= 1
832  floatdata(kdat+1)= -sngl(dps(m))
833  kdat=kdat+numcurv
834  DO l=1,offsetdimension
835  idx=kdat+l
836  intdata(idx)= ipar0+l
837  floatdata(idx)= sngl(wpp(m,l+l0))
838  idx=idx+offsetdimension
839  intdata(idx)= ipar0+l+offsetdimension
840  floatdata(idx)= sngl(wpn(m,l+l0))
841  END DO
842  kdat=kdat+2*offsetdimension
843  END IF
844  DO k=kdat0+1,kdat
845  IF (floatdata(k) /= 0.0) THEN ! copy non zero derivatives
846  integdatasize=integdatasize+1
847  intdata(integdatasize)=intdata(k)
848  floatdata(integdatasize)=floatdata(k)
849  ipmn=min(ipmn,intdata(integdatasize))
850  END IF
851  END DO
852  intdata(ldat0+1)=integdatasize-ldat0
853  ! check for local derivatives
854  ilcl=points(jpnt)%indLocal(m)
855  IF (ilcl <= 0) cycle
856  nlcl=intdata(ilcl)
857  IF (integdatasize+nlcl > maxdatasize) cycle
858  DO k=1,nlcl
859  IF (intdata(ilcl+k) <= numaddlocpar) THEN
860  integdatasize=integdatasize+1
861  intdata(integdatasize)=intdata(ilcl+k)+numcurv
862  floatdata(integdatasize)=floatdata(ilcl+k)
863  ipmn=min(ipmn,intdata(integdatasize))
864  END IF
865  END DO
866  intdata(ldat0+1)=integdatasize-ldat0
867  END DO
868 
869  END IF
870 
871  END DO
872 
873  IF (extseedpoint /= 0) THEN
874  IF (extseedmat(1) > 0.0d0) ipmn=1 ! external 'seed' for curvature
875 
876  nb=numcurv+numaddlocpar
877  npsd=nb+2*offsetdimension ! fit parameters
878  nplc=5+numaddlocpar ! local track parameters
879 
880  CALL gbljac(extseedpoint,1,ioff1) ! get transposed jacobian broken lines -> local
881  CALL devrot(nplc,diag,eigen,extseedmat,work,iwork)
882  forall(i=1:nplc,j=1:npsd) seed((j-1)*nplc+i)=sum(eigen((i-1)*nplc+1:i*nplc) &
883  *jactracktofit((j-1)*nplc+1:j*nplc)) ! eigen^T * jacTrackToFit
884  ip0=(ioff1-1)*offsetdimension
885 
886  DO i=1, nplc
887  IF (diag(i) > 0.0) THEN
888  IF (integdatasize+3+npsd > maxdatasize) return
889  ldat0=integdatasize
890  numdatablocks=numdatablocks+1 ! 'virtual' measurement
891  intdata(integdatasize+1)=3 ! length of item
892  intdata(integdatasize+2)=extseedpoint ! point
893  intdata(integdatasize+3)=i ! measurement
894  floatdata(integdatasize+1)=0.0 ! value
895  floatdata(integdatasize+2)=sngl(diag(i)) ! precision (1/error**2)
896  floatdata(integdatasize+3)=1.0 ! (additional) weight
897  integdatasize=integdatasize+3
898  DO j=1,npsd
899  IF (seed((j-1)*nplc+i) /= 0.0) THEN
900  integdatasize=integdatasize+1
901  IF (j > nb) THEN
902  intdata(integdatasize)=j+ip0
903  ELSE
904  intdata(integdatasize)=j
905  END IF
906  floatdata(integdatasize)=sngl(seed((j-1)*nplc+i))
907  END IF
908  END DO
909  intdata(ldat0+1)=integdatasize-ldat0
910  END IF
911  END DO
912  END IF
913 
914  ! check minimum parameter with non zero derivatives
915  IF (ipmn > numcurv) THEN ! curvature undefined
916  ipoff=numcurv
917  numcurv=0
918  ! correct parameter indices
919  IF (ipoff > 0) THEN
920  numfitpar=numfitpar-ipoff
921  integdatasize=0
922  DO i=1,numdatablocks
923  ltem=intdata(integdatasize+1)
924  forall(j=integdatasize+4:integdatasize+ltem) intdata(j)=intdata(j)-ipoff
925  integdatasize=integdatasize+ltem
926  END DO
927 
928  END IF
929  END IF
930  return
931  END SUBROUTINE gbldat
932 
940 
941  SUBROUTINE gblder(ipoint,idir,w,wj,dw)
942  IMPLICIT NONE
943 
944  INTEGER, INTENT(IN OUT) :: ipoint
945  INTEGER, INTENT(IN) :: idir
946  DOUBLE PRECISION, INTENT(OUT) :: w(2,2)
947  DOUBLE PRECISION, INTENT(OUT) :: wj(2,2)
948  DOUBLE PRECISION, INTENT(OUT) :: dw(2)
949 
950  DOUBLE PRECISION :: jm(2,2)
951  DOUBLE PRECISION :: sm(2,2)
952  DOUBLE PRECISION :: dv(2)
953  DOUBLE PRECISION :: det
954  INTEGER :: i
955  INTEGER :: j
956 
957  ! (parts of) jacobians to previous/next offset
958  IF (idir < 0) THEN
959  jm=points(ipoint)%prevJ
960  sm=-points(ipoint)%prevS
961  dv=points(ipoint)%prevd
962  ELSE
963  jm=points(ipoint)%nextJ
964  sm=points(ipoint)%nextS
965  dv=points(ipoint)%nextd
966  ENDIF
967 
968  det=sm(1,1)*sm(2,2)-sm(1,2)*sm(2,1) ! W
969  w(1,1)= sm(2,2)/det;w(1,2)=-sm(1,2)/det
970  w(2,1)=-sm(2,1)/det;w(2,2)= sm(1,1)/det
971  forall(i=1:2,j=1:2) wj(i,j)=sum(w(i,:)*jm(:,j)) ! W*J
972  forall(i=1:2) dw(i)=sum(dv(:)*w(i,:)) ! W*d
973 
974  return
975  END SUBROUTINE gblder
976 
978 
979  SUBROUTINE gblmat
980  IMPLICIT NONE
981  INTEGER :: i
982  INTEGER :: ij
983  INTEGER :: ij0
984  INTEGER :: ijsym
985  INTEGER :: ik
986  INTEGER :: j
987  INTEGER :: jk
988  INTEGER :: k
989  INTEGER :: ltem
990  INTEGER :: mpar2
991  INTEGER :: nbdr
992  INTEGER :: npar2
993 
994  DOUBLE PRECISION :: dval
995  DOUBLE PRECISION :: dwgh
996 
997  ! index in symmetric matrix
998  ijsym(i,j)=(i*i-i)/2+j
999 
1000  vecb(1:numfitpar)=0.0d0
1001  ! 'smart' clear
1002  IF (numfitpar > maxusedfitpar) THEN
1003  mpar2=(maxusedfitpar*maxusedfitpar+maxusedfitpar)/2
1004  mata(mpar2+1:npar2)=0.0d0
1005  maxusedfitpar=numfitpar
1006  END IF
1007 
1008  IF (lastnumfitpar > 0) THEN
1009  DO i=1,lastnumfitpar
1010  ij0=(i*i-i)/2
1011  forall(j=1:min(lastbordersize,i)) mata(ij0+j)=0.0d0 ! clear border
1012  forall(j=1:max(i-lastbandwidth,1):i) mata(ij0+j)=0.0d0 ! clear band
1013  END DO
1014  END IF
1015 
1016  integdatasize=0
1017  nbdr=numcurv+numaddlocpar
1018  DO i=1,numdatablocks
1019  ltem=intdata(integdatasize+1)
1020  dval=dble(floatdata(integdatasize+1)) ! value
1021  dwgh=dble(floatdata(integdatasize+3)*floatdata(integdatasize+2)) ! (total) weight
1022  DO j=integdatasize+4,integdatasize+ltem ! update matrix
1023  ij=intdata(j)
1024  vecb(ij)=vecb(ij)+dble(floatdata(j))*dval*dwgh
1025  DO k=integdatasize+4,j
1026  ik=intdata(k)
1027  jk=ijsym(ij,ik) ! parameter labels must be sorted: IK<=IJ !
1028  mata(jk)=mata(jk)+dble(floatdata(j))*dble(floatdata(k))*dwgh
1029  IF (ik > nbdr) bandsize=max(bandsize,ij-ik) ! update band width
1030  END DO
1031  END DO
1032  integdatasize=integdatasize+ltem
1033  END DO
1034 
1035  lastnumfitpar=numfitpar
1036  lastbordersize=nbdr
1037  lastbandwidth=bandsize
1038 
1039  return
1040  END SUBROUTINE gblmat
1041 
1047 
1048  SUBROUTINE gblch2(idwm,chi2,swgt)
1049  IMPLICIT NONE
1050  INTEGER :: i
1051  INTEGER :: j
1052  INTEGER :: ltem
1053  REAL :: brl
1054  REAL :: dwint
1055  REAL :: sres
1056  REAL :: val
1057  REAL :: wgh
1058 
1059  INTEGER, INTENT(IN) :: idwm
1060  REAL, INTENT(OUT) :: chi2
1061  REAL, INTENT(OUT) :: swgt
1062 
1063  dimension dwint(0:3) ! Integral(weight*normal_distribution)
1064  DATA dwint / 1.0, 0.8737, 0.9326, 0.8228 /
1065 
1066  chi2=0.0
1067  swgt=0.0
1068 
1069  integdatasize=0
1070  DO i=1,numdatablocks
1071  ltem=intdata(integdatasize+1)
1072  val=floatdata(integdatasize+1) ! value
1073  wgh=floatdata(integdatasize+3) ! down weighting
1074  brl=0.0 ! predction from broken lines
1075  DO j=integdatasize+4,integdatasize+ltem
1076  brl=brl+sngl(vecb(intdata(j))*floatdata(j))
1077  ENDDO
1078  sres=abs(brl-val)*sqrt(floatdata(integdatasize+2))
1079  chi2=chi2+sres*sres*wgh
1080  swgt=swgt+wgh
1081  ! down weighting
1082  IF (idwm == 1) THEN
1083  IF(sres < 4.6851) THEN ! Tukey
1084  wgh=(1.0-0.045558*sres*sres)**2 ! 1/4.6851**2
1085  ELSE
1086  wgh=0.0
1087  END IF
1088  ELSE IF (idwm == 2) THEN
1089  IF (sres < 1.345) THEN ! Huber
1090  wgh=1.0
1091  ELSE
1092  wgh=1.345/sres
1093  END IF
1094  ELSE IF (idwm == 3) THEN
1095  wgh=1.0/(1.0+(sres/2.3849)**2) ! Cauchy
1096  END IF
1097  floatdata(integdatasize+3)=wgh
1098  integdatasize=integdatasize+ltem
1099  END DO
1100 
1101  chi2=chi2/dwint(lastdownweightmethod) ! renormalize Chi2
1102  lastdownweightmethod=idwm
1103 
1104  return
1105  END SUBROUTINE gblch2
1106 
1110 
1111  SUBROUTINE gblmp2(iret)
1112  IMPLICIT NONE
1113  INTEGER :: i
1114  INTEGER :: igbl
1115  INTEGER :: im
1116  INTEGER :: ipnt
1117  INTEGER :: j
1118  INTEGER :: ltem
1119  INTEGER :: ngbl
1120  REAL :: sig
1121  REAL :: derlc
1122 
1123  INTEGER, INTENT(OUT) :: iret
1124 
1125  dimension derlc(maxfitpar)
1126 
1127  iret=0
1128  CALL gblprp
1129  IF (trajlevel < 3) return
1130 
1131  integdatasize=0
1132  DO i=1,numdatablocks
1133  ltem=intdata(integdatasize+1)
1134  ipnt=intdata(integdatasize+2)
1135  im =intdata(integdatasize+3)
1136  ! local derivatives
1137  derlc(1:numfitpar)=0.0
1138  forall(j=integdatasize+4:integdatasize+ltem) derlc(intdata(j))=floatdata(j)
1139  igbl=0
1140  ngbl=0
1141  IF (im > 0) igbl=points(ipnt)%indGlobal(im)
1142  IF (igbl > 0) ngbl=intdata(igbl)
1143  sig=1.0/sqrt(floatdata(integdatasize+2))
1144 
1145  CALL mille(numfitpar,derlc,ngbl,floatdata(igbl+1),intdata(igbl+1), &
1146  floatdata(integdatasize+1),sig) ! add data
1147 
1148  integdatasize=integdatasize+ltem
1149  END DO
1150 
1151  CALL endle ! complete, write record (many sets)
1152  iret=numdatablocks
1153 
1154  return
1155  END SUBROUTINE gblmp2
1156 
1164 
1165  SUBROUTINE gbljac(ipoint,itrans,ioff1)
1166  IMPLICIT NONE
1167  INTEGER :: i
1168  INTEGER :: idir
1169  INTEGER :: indx
1170  INTEGER :: ioff2
1171  INTEGER :: ip1
1172  INTEGER :: ip2
1173  INTEGER :: j
1174  INTEGER :: joff
1175  INTEGER :: jpnt
1176  INTEGER :: k
1177  INTEGER :: koff
1178  INTEGER :: l
1179  INTEGER :: l0
1180  INTEGER :: mp
1181  INTEGER :: np
1182 
1183  INTEGER, INTENT(IN) :: ipoint
1184  INTEGER, INTENT(IN) :: itrans
1185  INTEGER, INTENT(OUT) :: ioff1
1186 
1187  DOUBLE PRECISION :: w(2,2)
1188  DOUBLE PRECISION :: wj(2,2)
1189  DOUBLE PRECISION :: dw(2)
1190  DOUBLE PRECISION :: wp(2,2)
1191  DOUBLE PRECISION :: wjp(2,2)
1192  DOUBLE PRECISION :: wn(2,2)
1193  DOUBLE PRECISION :: wjn(2,2)
1194  DOUBLE PRECISION :: dwp(2)
1195  DOUBLE PRECISION :: dwn(2)
1196  DOUBLE PRECISION :: wjs(2,2)
1197  DOUBLE PRECISION :: wji(2,2)
1198  DOUBLE PRECISION :: wip(2,2)
1199  DOUBLE PRECISION :: win(2,2)
1200  DOUBLE PRECISION :: dip(2)
1201  DOUBLE PRECISION :: din(2)
1202  DOUBLE PRECISION :: wpp(2,2)
1203  DOUBLE PRECISION :: wpn(2,2)
1204  DOUBLE PRECISION :: dpp(2)
1205  DOUBLE PRECISION :: dpn(2)
1206  DOUBLE PRECISION :: det
1207  DOUBLE PRECISION :: sgn
1208 
1209  indx(i,j)=((i-1)*np+j)*(1-itrans)+((j-1)*mp+i)*itrans
1210 
1211  np=numcurv+numaddlocpar+2*offsetdimension
1212  mp=5+numaddlocpar
1213  jpnt=iabs(ipoint)
1214 
1215  joff=points(jpnt)%offset
1216 
1217  jactracktofit(1:mp*np)=0.0d0 ! reset Jacobi matrix
1218 
1219  l0=firstactivecoord
1220  IF (joff > 0) THEN
1221  ! at offset
1222  IF (ipoint > 0) THEN ! forward
1223  idir=2
1224  IF (joff == numoffsets) idir=-1
1225  ELSE ! backward
1226  idir=-1
1227  IF (joff == 1) idir=2
1228  END IF
1229  koff=joff+2*iabs(idir)-3 ! 2nd offset for slope measurement
1230  CALL gblder(jpnt,idir,w,wj,dw)
1231 
1232  IF (idir > 0) THEN
1233  ioff1=joff
1234  ioff2=koff
1235  ip1=numcurv+numaddlocpar
1236  ip2=numcurv+numaddlocpar+offsetdimension
1237  sgn=1.0
1238  ELSE
1239  ioff2=joff
1240  ioff1=koff
1241  ip2=numcurv+numaddlocpar
1242  ip1=numcurv+numaddlocpar+offsetdimension
1243  sgn=-1.0
1244  END IF
1245 
1246  DO l=1,offsetdimension
1247  jactracktofit(indx(2 ,ip1+l))=-sgn*wj(1,l+l0)
1248  jactracktofit(indx(3 ,ip1+l))=-sgn*wj(2,l+l0)
1249  jactracktofit(indx(2 ,ip2+l))= sgn*w(1,l+l0)
1250  jactracktofit(indx(3 ,ip2+l))= sgn*w(2,l+l0)
1251  jactracktofit(indx(3+l,ip1+l))= 1.0d0
1252  END DO
1253  IF (numcurv > 0) THEN ! correct for curvature
1254  jactracktofit(indx(2,1)) = -sgn*dw(1)
1255  jactracktofit(indx(3,1)) = -sgn*dw(2)
1256  END IF
1257 
1258  ELSE
1259  ! between offsets
1260  CALL gblder(jpnt,-1,wp,wjp,dwp)
1261  CALL gblder(jpnt, 2,wn,wjn,dwn)
1262  ioff1=-joff-1
1263  ioff2=-joff
1264  ip1=numcurv+numaddlocpar
1265  ip2=numcurv+numaddlocpar+offsetdimension
1266 
1267  wjs=wjp+wjn ! W- * J- + W+ * J+
1268  det=wjs(1,1)*wjs(2,2)-wjs(1,2)*wjs(2,1) ! (W- * J- + W+ * J+)^-1 (=N)
1269  wji(1,1)= wjs(2,2)/det;wji(1,2)=-wjs(1,2)/det
1270  wji(2,1)=-wjs(2,1)/det;wji(2,2)= wjs(1,1)/det
1271  ! derivatives for u_int
1272  forall(i=1:2,j=1:2) wip(i,j)=sum(wji(i,:)*wp(:,j)) ! N * W-
1273  forall(i=1:2,j=1:2) win(i,j)=sum(wji(i,:)*wn(:,j)) ! N * W+
1274  forall(i=1:2) dip(i)=sum(wji(i,:)*dwp(:)) ! N * W- * d-
1275  forall(i=1:2) din(i)=sum(wji(i,:)*dwn(:)) ! N * W+ * d+
1276  ! derivatives for alpha_int
1277  forall(i=1:2,j=1:2) wpp(i,j)=sum(wjn(i,:)*wip(:,j)) ! W+ * J+ * N * W-
1278  forall(i=1:2,j=1:2) wpn(i,j)=sum(wjp(i,:)*win(:,j)) ! W- * J- * N * W+
1279  forall(i=1:2) dpp(i)=sum(wjn(i,:)*dip(:)) ! W+ * J+ * N * W- * d-
1280  forall(i=1:2) dpn(i)=sum(wjp(i,:)*din(:)) ! W- * J- * N * W+ * d+
1281 
1282  DO l=1,offsetdimension
1283  ! du_int/du-
1284  jactracktofit(indx(4,ip1+l))= wip(1,l+l0)
1285  jactracktofit(indx(5,ip1+l))= wip(2,l+l0)
1286  ! du_int/du+
1287  jactracktofit(indx(4,ip2+l))= win(1,l+l0)
1288  jactracktofit(indx(5,ip2+l))= win(2,l+l0)
1289  ! dalpha_int/du-
1290  jactracktofit(indx(2,ip1+l))=-wpp(1,l+l0)
1291  jactracktofit(indx(3,ip1+l))=-wpp(2,l+l0)
1292  ! dalpha_int/du+
1293  jactracktofit(indx(2,ip2+l))= wpn(1,l+l0)
1294  jactracktofit(indx(3,ip2+l))= wpn(2,l+l0)
1295  END DO
1296  IF (numcurv > 0) THEN ! correct for curvature
1297  ! du_int/dQbyP
1298  jactracktofit(indx(4,1)) =-dip(1)-din(1)
1299  jactracktofit(indx(5,1)) =-dip(2)-din(2)
1300  ! dalpha_int/dQbyP
1301  jactracktofit(indx(2,1)) = dpp(1)-dpn(1)
1302  jactracktofit(indx(3,1)) = dpp(2)-dpn(2)
1303  END IF
1304 
1305  END IF
1306 
1307  IF (numcurv > 0) THEN ! curvature
1308  jactracktofit(indx(1,1))=1.0d0
1309  END IF
1310 
1311  DO k=1, numaddlocpar
1312  jactracktofit(indx(k+5,numcurv+k))= 1.0d0 ! local parameters
1313  END DO
1314 
1315  return
1316  END SUBROUTINE gbljac
1317 
1327 
1328  SUBROUTINE gblres(ipoint,dpar,dcov)
1329  IMPLICIT NONE
1330  INTEGER :: i
1331  INTEGER :: io
1332  INTEGER :: ioff1
1333  INTEGER :: ip0
1334  INTEGER :: j
1335  INTEGER :: jpnt
1336  INTEGER :: k
1337  INTEGER :: ko
1338  INTEGER :: mp
1339  INTEGER :: mp2
1340  INTEGER :: nb
1341  INTEGER :: nb2
1342  INTEGER :: np
1343 
1344  INTEGER, INTENT(IN OUT) :: ipoint
1345  DOUBLE PRECISION, INTENT(OUT) :: dpar(*)
1346  DOUBLE PRECISION, INTENT(OUT) :: dcov(*)
1347 
1348  DOUBLE PRECISION :: caux(maxtracksymmatsize)
1349  DOUBLE PRECISION :: baux(maxtrackpar)
1350 
1351  nb=numcurv+numaddlocpar
1352  nb2=(nb*nb+nb)/2
1353  np=nb+2*offsetdimension
1354  mp=5+numaddlocpar
1355  mp2=(mp*mp+mp)/2
1356 
1357  jpnt=iabs(ipoint)
1358  IF (jpnt < 1.OR.jpnt > numpoints) THEN
1359  IF (printlevel > 0) print *, ' GBLRES invalid point ', ipoint
1360  dpar(1:mp)=0.0d0
1361  dcov(1:mp2)=0.0d0
1362  return
1363  END IF
1364 
1365  IF (trajlevel < 4) return ! fit not yet performed or failed
1366 
1367  CALL gbljac(ipoint,0,ioff1) ! get jacobian broken lines -> local
1368  ! get compressed covariance matrix, result vector
1369  baux(1:nb)=vecb(1:nb) ! border
1370  caux(1:nb2)=mata(1:nb2) ! border
1371  ip0=(ioff1-1)*offsetdimension
1372  k=nb2
1373  DO i=nb+1,np
1374  io=i+ip0
1375  ko=(io*io-io)/2
1376  baux(i)=vecb(io) ! band part
1377  DO j=1,nb ! mixed part
1378  k=k+1
1379  caux(k)=mata(ko+j)
1380  END DO
1381  ko=ko+ip0
1382  DO j=nb+1,i ! band part
1383  k=k+1
1384  caux(k)=mata(ko+j)
1385  END DO
1386  END DO
1387 
1388  CALL dbgax(jactracktofit,baux,dpar,mp,np)
1389  CALL dbavat(caux,jactracktofit,dcov,np,mp)
1390 
1391  return
1392  END SUBROUTINE gblres
1393 
1403 
1404  SUBROUTINE gbladx(ipoint,dprc)
1405  IMPLICIT NONE
1406  INTEGER :: jpnt
1407  INTEGER :: mp
1408  INTEGER :: mp2
1409 
1410  INTEGER, INTENT(IN) :: ipoint
1411  DOUBLE PRECISION, INTENT(IN) :: dprc(*)
1412 
1413 
1414  jpnt=iabs(ipoint)
1415  IF (jpnt < 1.OR.jpnt > numpoints) THEN
1416  IF (printlevel > 0) print *, ' GBLADX invalid point ', ipoint
1417  return
1418  END IF
1419 
1420  IF (trajlevel >= 3) return ! fit already prepared or performed
1421 
1422  mp=5+numaddlocpar
1423  mp2=(mp*mp+mp)/2
1424  extseedpoint=ipoint
1425  extseedmat(1:mp2)=dprc(1:mp2)
1426 
1427  return
1428  END SUBROUTINE gbladx
1429 
1435  SUBROUTINE gblinu(a,b)
1436  IMPLICIT NONE
1437 
1438  DOUBLE PRECISION, INTENT(IN) :: a(5,5)
1439  DOUBLE PRECISION, INTENT(OUT) :: b(5,5)
1440  DOUBLE PRECISION :: a33(3,3)
1441  DOUBLE PRECISION :: a23(2,3)
1442  DOUBLE PRECISION :: a22(2,2)
1443  DOUBLE PRECISION :: det
1444 
1445  ! A33 = A(1..3,1..3)^-1 * det(A(1..3,1..3)
1446  a33(1,1) = a(2,2)*a(3,3)-a(2,3)*a(3,2)
1447  a33(2,1) = a(2,3)*a(3,1)-a(2,1)*a(3,3)
1448  a33(3,1) = a(2,1)*a(3,2)-a(2,2)*a(3,1)
1449  a33(1,2) = a(1,3)*a(3,2)-a(1,2)*a(3,3)
1450  a33(2,2) = a(1,1)*a(3,3)-a(1,3)*a(3,1)
1451  a33(3,2) = a(1,2)*a(3,1)-a(1,1)*a(3,2)
1452  a33(1,3) = a(1,2)*a(2,3)-a(1,3)*a(2,2)
1453  a33(2,3) = a(1,3)*a(2,1)-a(1,1)*a(2,3)
1454  a33(3,3) = a(1,1)*a(2,2)-a(1,2)*a(2,1)
1455  det=a(1,1)*a33(1,1)+a(1,2)*a33(2,1)+a(1,3)*a33(3,1)
1456  ! A23 = A(4..5,1..3) * A33 / det(A(1..3,1..3))
1457  a23(1,1) = (a(4,1)*a33(1,1)+a(4,2)*a33(2,1)+a(4,3)*a33(3,1))/det
1458  a23(1,2) = (a(4,1)*a33(1,2)+a(4,2)*a33(2,2)+a(4,3)*a33(3,2))/det
1459  a23(1,3) = (a(4,1)*a33(1,3)+a(4,2)*a33(2,3)+a(4,3)*a33(3,3))/det
1460  a23(2,1) = (a(5,1)*a33(1,1)+a(5,2)*a33(2,1)+a(5,3)*a33(3,1))/det
1461  a23(2,2) = (a(5,1)*a33(1,2)+a(5,2)*a33(2,2)+a(5,3)*a33(3,2))/det
1462  a23(2,3) = (a(5,1)*a33(1,3)+a(5,2)*a33(2,3)+a(5,3)*a33(3,3))/det
1463  ! A22 = A(4..5,4..5) - A23 * A((1..3,4..5)
1464  a22(1,1) = a(4,4)-a23(1,1)*a(1,4)-a23(1,2)*a(2,4)-a23(1,3)*a(3,4)
1465  a22(1,2) = a(4,5)-a23(1,1)*a(1,5)-a23(1,2)*a(2,5)-a23(1,3)*a(3,5)
1466  a22(2,1) = a(5,4)-a23(2,1)*a(1,4)-a23(2,2)*a(2,4)-a23(2,3)*a(3,4)
1467  a22(2,2) = a(5,5)-a23(2,1)*a(1,5)-a23(2,2)*a(2,5)-a23(2,3)*a(3,5)
1468  ! (4..5,4..5) of A^-1 = A22^-1
1469  det = a22(1,1)*a22(2,2)-a22(1,2)*a22(2,1)
1470  b(4,4) = a22(2,2)/det
1471  b(4,5) =-a22(1,2)/det
1472  b(5,4) =-a22(2,1)/det
1473  b(5,5) = a22(1,1)/det
1474  ! (4..5,1..3) of A^-1 = -A22^-1 * A23
1475  b(4,1) = -b(4,4)*a23(1,1)-b(4,5)*a23(2,1)
1476  b(4,2) = -b(4,4)*a23(1,2)-b(4,5)*a23(2,2)
1477  b(4,3) = -b(4,4)*a23(1,3)-b(4,5)*a23(2,3)
1478  b(5,1) = -b(5,4)*a23(1,1)-b(5,5)*a23(2,1)
1479  b(5,2) = -b(5,4)*a23(1,2)-b(5,5)*a23(2,2)
1480  b(5,3) = -b(5,4)*a23(1,3)-b(5,5)*a23(2,3)
1481  return
1482  END SUBROUTINE gblinu
1483 END MODULE gbltraj
1484 
Generated on Mon Jan 21 2013 07:17:06 for GeneralBrokenLines by   doxygen 1.8.1