10Â¥: Originally posted by wlfc at 2021-08-11 16:09:33
ʹÓà 3Â¥¡¢4Â¥ µÄ¹«Ê½½øÐÐÇó½â¡£

ÓÐÁ½¸öÄѵ㣺
1¡¢Î¢·Ö·½³ÌÇó½âʱҪ´«µÝÖмä±äÁ¿lamb, pa¡£
2¡¢lambvΪδ֪Öмä±äÁ¿£¬ÐèÒªÇó½â¡£
µÚÒ»¸öÄѵã¶ÔLu½Å±¾À´Ëµ²»ÊÇÎÊÌ⣬gslÖ§³ÖµÄLuÀ©Õ¹Êýѧ¿âÖÐÇó½â΢·Ö·½³ÌµÄº¯Êý ...

7Â¥: Originally posted by wlfc at 2021-08-09 08:09:41
³¢ÊÔÓÃOpenLu£¨¿É´Ówww.forcal.netÏÂÔØ£©Çó½â£¬Óë´ó¼Ò̽ÌÖ¡£
ʹÓà 1Â¥ ¸ø³öµÄµÚÒ»ÖÖ¹«Ê½¡£Î´ÖªÖмä±äÁ¿lambvÓýⷽ³Ì·¨Çó½â£¬µ«²»ÖªÕâÖÖ·½·¨ÓÐûÓÐÎÊÌ⣿
ÈÔÀàËÆÓÚ 4Â¥ ʹÓò¿·ÖÊý¾ÝÄâºÏ£¬ÒÔ¼õÉÙÔËÐÐʱ¼ä¡£

Lu ...

!!!using["luopt","math"]; //ʹÓÃÃüÃû¿Õ¼ä

f(  t,lambv,dlambv,i
    : lamb
    : t_A, t_lamb, fac1,fac2,fac3,fac4,fac5,fac6,fac7) =
    lamb=t_lamb[i]+(t_lamb[i+1]-t_lamb[i])*[(t-t_A[i])/(t_A[i+1]-t_A[i])],
    dlambv=((1.0/3/fac7)*(fac1*((lamb/lambv)^fac2-(lambv/lamb)^(0.5*fac2))+fac3*((lamb/lambv)^fac4-(lambv/lamb)^(0.5*fac4))+fac5*((lamb/lambv)^fac6-(lambv/lamb)^(0.5*fac6)))),
    0;  //³É¹¦·µ»Ø0

Ä¿±êº¯Êý(_fac1,_fac2,_fac3,_fac4,_fac5,_fac6,_fac7,lambv0
    : i,s,tf,pa,lamb,lambv
    : t_A, t_lamb, t_pa, t_p, max, fac1,fac2,fac3,fac4,fac5,fac6,fac7)=
{
    fac1=_fac1, fac2=_fac2, fac3=_fac3, fac4=_fac4, fac5=_fac5, fac6=_fac6, fac7=_fac7,  //´«µÝÓÅ»¯±äÁ¿
    //×îºóÒ»¸ö²ÎÊý50±íʾgsl_odeº¯ÊýÔÚ¼ÆËãʱ,×î¶àÑ­»·¼ÆËã50´Î,ÕâÑù¿ÉÒÔÌá¸ßËÙ¶È
    tf=gsl_ode[@f, nil, nil, t_A, ra1(lambv0),  1e-6, 1e-6, gsl_rkf45, 1e-6,50],
    i=-1, s=0, while{++i<max,
        pa=t_pa(i), lamb=t_lamb(i), lambv=tf(i,1),
        s=s+[(pa+fac1*((lamb/lambv)^fac2/lamb-(lambv/lamb)^(fac2*0.5)/lamb)+fac3*((lamb/lambv)^fac4/lamb-(lambv/lamb)^(fac4*0.5)/lamb)+fac5*((lamb/lambv)^fac6/lamb-(lambv/lamb)^(fac6*0.5)/lamb)) - t_p(i)]^2.0
    },
    s
};

main(: tArray : t_A, t_lamb, t_pa, t_p, max)=
{
    tArray=matrix{ //´æ·ÅʵÑéÊý¾Ý//t,lamb,pa,p
        "0        0.959926667        -0.218715994        -0.060966222
0.25        0.963063333        -0.200840559        -0.046836278
0.5        0.96612        -0.183553472        -0.035042256
0.75        0.969016667        -0.167289626        -0.024438988
1        0.97171        -0.152268899        -0.01706765
1.25        0.97409        -0.139075651        -0.010939712
1.5        0.976123333        -0.127862697        -0.006795087
1.75        0.977786667        -0.118729759        -0.002919842
2        0.979000001        -0.112089924        -0.001981819
2.25        0.979743333        -0.108031322        -0.001163676
2.5        0.98001        -0.106577018        -0.003491456
2.75        0.979766667        -0.107904033        -0.009723697"
    },
    len[tArray,0,&max], t_A=tArray(all:0).reshape(), t_lamb=tArray(all:1).reshape(), t_pa=tArray(all:2).reshape(), t_p=tArray(all:3).reshape(), //ÓÃlenº¯ÊýÈ¡¾ØÕóµÄÐÐÊý£¬t_AµÈÈ¡¾ØÕóµÄÁÐ
    Opt1[@Ä¿±êº¯Êý,  optwaysimdeep, optwayconfra] //Opt1º¯ÊýÈ«¾ÖÓÅ»¯
};

!!!using["luopt","math","win"]; //ʹÓÃÃüÃû¿Õ¼ä

f(  t,lambv,dlambv,i
    : lamb
    : t_A, t_lamb, fac1,fac2,fac3,fac4,fac5,fac6,fac7) =
    lamb=t_lamb[i]+(t_lamb[i+1]-t_lamb[i])*[(t-t_A[i])/(t_A[i+1]-t_A[i])],
    dlambv=((1.0/3/fac7)*(fac1*((lamb/lambv)^fac2-(lambv/lamb)^(0.5*fac2))+fac3*((lamb/lambv)^fac4-(lambv/lamb)^(0.5*fac4))+fac5*((lamb/lambv)^fac6-(lambv/lamb)^(0.5*fac6)))),
    0;  //³É¹¦·µ»Ø0

get(t_pp,t_lambv, _fac1,_fac2,_fac3,_fac4,_fac5,_fac6,_fac7,lambv0
    : i,tf,pa,lamb,lambv
    : t_A, t_lamb, t_pa, t_p, max, fac1,fac2,fac3,fac4,fac5,fac6,fac7)=
{
    fac1=_fac1, fac2=_fac2, fac3=_fac3, fac4=_fac4, fac5=_fac5, fac6=_fac6, fac7=_fac7,  //´«µÝÓÅ»¯±äÁ¿
    //×îºóÒ»¸ö²ÎÊý50±íʾgsl_odeº¯ÊýÔÚ¼ÆËãʱ,×î¶àÑ­»·¼ÆËã50´Î,ÕâÑù¿ÉÒÔÌá¸ßËÙ¶È
    tf=gsl_ode[@f, nil, nil, t_A, ra1(lambv0),  1e-6, 1e-6, gsl_rkf45, 1e-6,50],
    t_lambv.copy[tf(all,1)],
    i=-1, while{++i<max,
        pa=t_pa(i), lamb=t_lamb(i), lambv=t_lambv(i),
        t_pp[i]=(pa+fac1*((lamb/lambv)^fac2/lamb-(lambv/lamb)^(fac2*0.5)/lamb)+fac3*((lamb/lambv)^fac4/lamb-(lambv/lamb)^(fac4*0.5)/lamb)+fac5*((lamb/lambv)^fac6/lamb-(lambv/lamb)^(fac6*0.5)/lamb))
    }
};

init( main
    :  tArray, tf, i, k, t_pp,t_lambv
    : t_A, t_lamb, t_pa, t_p, max)=
{
    tArray=matrix{ //´æ·ÅʵÑéÊý¾Ý//t,lamb,pa,p
        "0        0.959926667        -0.218715994        -0.060966222
0.25        0.963063333        -0.200840559        -0.046836278
0.5        0.96612        -0.183553472        -0.035042256
0.75        0.969016667        -0.167289626        -0.024438988
1        0.97171        -0.152268899        -0.01706765
1.25        0.97409        -0.139075651        -0.010939712
1.5        0.976123333        -0.127862697        -0.006795087
1.75        0.977786667        -0.118729759        -0.002919842
2        0.979000001        -0.112089924        -0.001981819
2.25        0.979743333        -0.108031322        -0.001163676
2.5        0.98001        -0.106577018        -0.003491456
2.75        0.979766667        -0.107904033        -0.009723697"
    },
    len[tArray,0,&max], t_A=tArray(all:0).reshape(), t_lamb=tArray(all:1).reshape(), t_pa=tArray(all:2).reshape(), t_p=tArray(all:3).reshape(), //ÓÃlenº¯ÊýÈ¡¾ØÕóµÄÐÐÊý£¬t_AµÈÈ¡¾ØÕóµÄÁÐ
    t_pp=new[real_s,max], t_lambv=new[real_s,max],

    get[t_pp, t_lambv, 7.148976005697928e-003  ,  -14.5708961980305     ,    -1.059296609366114    ,    0.2342635320112431    ,    -9.68049051324744e-005  ,  -31.79361445870142     ,   -2.026215042289879     ,   1.184269268536115  ],

    cwAttach[typeSplit], cwResizePlots(1,2,2), //×ó¶þÓÒ¶þ·ÖÁÑ
    k=cwAddCurve{t_A, t_p, max, 0},    //¸ø0×ÓͼÌí¼ÓÇúÏß
    cwSetScatter(k,0),                 //ÉèÖûæÖƵã
    cwSetDataLineSize(5, k, 0),        //ÉèÖõãµÄ´óС
    cwAddCurve{t_A, t_pp, max, 0},     //¸ø0×ÓͼÌí¼ÓÇúÏß
    cwAddCurve{t_A, t_lamb, max, 1},   //¸ø1×ÓͼÌí¼ÓÇúÏß
    cwAddCurve{t_A, t_pa, max, 2},     //¸ø2×ÓͼÌí¼ÓÇúÏß
    cwAddCurve{t_A, t_lambv, max, 3}   //¸ø3×ÓͼÌí¼ÓÇúÏß
};
ChartWnd[@init];