45 #ifndef CLIPPER_RESOL_TARGETFN
46 #define CLIPPER_RESOL_TARGETFN
48 #include "resol_basisfn.h"
49 #include "hkl_datatypes.h"
239 ftype omeg = (omegaa>0.05) ? omegaa : (0.05*exp(omegaa/0.05-1.0));
240 return 0.5 * ( sqrt( 4.0*omeg*omeg + 1.0 ) - 1.0 ) / omeg;
284 hkl_data = &hkl_data_;
291 if ( !data[ih].missing() ) {
298 result.
r = result.
dr = result.
dr2 = 0.0;
308 hkl_data = &hkl_data_;
315 if ( !data[ih].missing() ) {
316 ftype d = fh - pow(
ftype(data[ih].E()), power );
321 result.
r = result.
dr = result.
dr2 = 0.0;
330 hkl_data1 = &hkl_data1_;
331 hkl_data2 = &hkl_data2_;
337 const T1& ft1 = (*hkl_data1)[ih];
338 const T2& ft2 = (*hkl_data2)[ih];
339 if ( !ft1.missing() && !ft2.missing() ) {
341 const ftype f1 = pow( ft1.f(), 2 ) / eps;
342 const ftype f2 = pow( ft2.f(), 2 ) / eps;
343 const ftype d = fh*f1 - f2;
344 result.
r = d * d / f1;
346 result.
dr2 = 2.0 * f1;
348 result.
r = result.
dr = result.
dr2 = 0.0;
357 hkl_data1 = &hkl_data1_;
358 hkl_data2 = &hkl_data2_;
364 result.
r = result.
dr = result.
dr2 = 0.0;
365 const T1& ft1 = (*hkl_data1)[ih];
366 const T2& ft2 = (*hkl_data2)[ih];
367 if ( !ft1.missing() && !ft2.missing() )
368 if ( ft1.f() > 1.0e-6 && ft2.f() > 1.0e-6 ) {
370 const ftype f1 = pow( ft1.f(), 2 ) / eps;
371 const ftype f2 = pow( ft2.f(), 2 ) / eps;
373 const ftype d = fh + log(f1) - log(f2);
374 result.
r = w * d * d;
375 result.
dr = 2.0 * w * d;
376 result.
dr2 = 2.0 * w;
385 hkl_data1 = &hkl_data1_;
386 hkl_data2 = &hkl_data2_;
392 const T1& ft1 = (*hkl_data1)[ih];
393 const T2& ft2 = (*hkl_data2)[ih];
394 if ( !ft1.missing() && !ft2.missing() ) {
396 const ftype f1 = ft1.I() / eps;
397 const ftype f2 = ft2.I() / eps;
398 const ftype d = fh*f1 - f2;
399 result.
r = d * d / f1;
401 result.
dr2 = 2.0 * f1;
403 result.
r = result.
dr = result.
dr2 = 0.0;
412 hkl_data1 = &hkl_data1_;
413 hkl_data2 = &hkl_data2_;
419 result.
r = result.
dr = result.
dr2 = 0.0;
420 const T1& ft1 = (*hkl_data1)[ih];
421 const T2& ft2 = (*hkl_data2)[ih];
422 if ( !ft1.missing() && !ft2.missing() )
423 if ( ft1.I() > 1.0e-6 && ft2.I() > 1.0e-6 ) {
425 const ftype f1 = ft1.I() / eps;
426 const ftype f2 = ft2.I() / eps;
428 const ftype d = fh + log(f1) - log(f2);
429 result.
r = w * d * d;
430 result.
dr = 2.0 * w * d;
431 result.
dr2 = 2.0 * w;
440 hkl_data = &hkl_data_;
447 const ftype two(2.0);
448 if ( !data[ih].missing() ) {
450 ftype d = fsq * fh - 1.0;
451 result.
r = d * d / fsq;
453 result.
dr2 = two * fsq;
455 result.
r = result.
dr = result.
dr2 = 0.0;
475 if ( eodata[ih].missing() || ecdata[ih].missing() ) {
476 result.
r = result.
dr = result.
dr2 = 0.0;
478 ftype eo = eodata[ih].E();
479 ftype ec = ecdata[ih].E();
480 ftype omeg = (omegaa>0.05) ? omegaa : (0.05*exp(omegaa/0.05-1.0));
481 ftype sigmaa = 0.5*(sqrt(4*omeg*omeg+1)-1)/omeg;
482 ftype dx = 2.0 * eo * ec;
484 ftype t0 = 1.0/(1-sigmaa*sigmaa) + 0.5*log((1-sigmaa*sigmaa));
486 ftype t2 = pow(1-sigmaa*sigmaa,2)/(1+sigmaa*sigmaa);
488 result.
r = 1.0*t0 - log( cosh( x/2 ) );
489 result.
dr = 1.0*t1 - dx*0.5*tanh( x/2 );
490 result.
dr2 = 1.0*t2 - dx*dx*0.25*(1.0 - pow(tanh(x/2),2) );
496 if ( omegaa < 0.05 ) {
497 ftype dy = exp( omegaa/0.05 ) / exp( 1.0 );
498 ftype dy2 = exp( omegaa/0.05 ) / ( 0.05*exp( 1.0 ) );
499 result.
dr2 = result.
dr*dy2 + result.
dr2*dy*dy;
500 result.
dr = result.
dr*dy;
521 if ( eodata[ih].missing() || ecdata[ih].missing() ) {
522 result.
r = result.
dr = result.
dr2 = 0.0;
524 ftype eo = eodata[ih].E();
525 ftype ec = ecdata[ih].E();
526 ftype sigmaa = (sigmaa0>0.99)?(0.99):((sigmaa0<0.01)?0.01:sigmaa0);
527 ftype dx = 2.0 * eo * ec;
528 ftype x = dx * sigmaa/(1-sigmaa*sigmaa);
529 ftype t0 = 1.0/(1-sigmaa*sigmaa) + 0.5*log((1-sigmaa*sigmaa));
531 ftype t2 = pow(1-sigmaa*sigmaa,2)/(1+sigmaa*sigmaa);
533 result.
r = 1.0*t0 - log( cosh( x/2 ) );
534 result.
dr = 1.0*t1 - dx*0.5*tanh( x/2 );
535 result.
dr2 = 1.0*t2 - dx*dx*0.25*(1.0 - pow(tanh(x/2),2) );
541 ftype ds = (1+sigmaa*sigmaa)/pow(1-sigmaa*sigmaa,2);
542 ftype ds2 = 2*sigmaa*(3+sigmaa*sigmaa)/pow(1-sigmaa*sigmaa,3);
543 result.
dr2 = result.
dr*ds2 + result.
dr2*ds*ds;
544 result.
dr = result.
dr*ds;