; ______________________________________________________________;
;                     MODEL10.PRO
; -Bootstrap example, where minimum and total amplitude solved
; numerically.
; ______________________________________________________________;
PRO FU,X,A,F,PDER
F=0.+A(0)+A(1)*COS(X)+A(2)*SIN(X)  & PDER=DBLARR(N_ELEMENTS(X),3)
PDER(*,0)=1.D0  & PDER(*,1)=COS(X) & PDER(*,2)=SIN(X)
RETURN & END    
; ______________________________________________________________;
PRO MITTAA,XX,GG,AMP,TMIN
; - Solves the total amplitude and first minimum
J=WHERE((GG LT SHIFT(GG,1)) AND (GG LT SHIFT(GG,-1))$
AND     (XX NE MIN(XX))     AND (XX NE MAX(XX))    ) 
TMIN=MIN(XX(J(0)))                 & AMP=MAX(GG)-MIN(GG)
RETURN & END
; ______________________________________________________________;
; -Simulated data
N =10.D0 & READ,'Give ..... n = ',N
Z =10.D0 & READ,'Give A/sigma = ',Z
M =10.D0 & READ,'Give ..... M = ',M 
B_1=1.D0 & READ,'Give ... B_1 = ',B_1
C_1=1.D0 & READ,'Give ... C_1 = ',C_1
T =2.D0*RANDOMU(S,N)         & T =2.D0*!PI*T
Y =M+B_1*COS(T)+C_1*SIN(T)   & A =MAX(Y)-MIN(Y)
E =(A/Z)*RANDOMN(S,N)        & Y=Y+E   & E =(A/Z)*RANDOMN(S,N)
J =SORT(T) & T=T(J) & Y=Y(J) & E=E(J) ; Sorted in time
; ______________________________________________________________;
!P.NOERASE=1  & !X.TICKS=2. & !Y.TICKS=2. & !P.SYMSIZE=1.5
!X.STYLE=1    & !Y.STYLE=1  & !P.CHARSIZE=1.5
; ______________________________________________________________; 
ROUNDS=10.D0 & READ,'Bootstrap samples     = ',ROUNDS
AIKAA =5.D0  & READ,'Seconds between plots = ',AIKAA 
W=E^(-2.D0)  & BETA=DBLARR(3)
G=CURVEFIT(T,Y,W,BETA,SIGMAA,FUNCTION_NAME='FU') & EPSILON=Y-G
; ______________________________________________________________; 
X1=MIN(T) & X2=MAX(T) & Y1=0.9*MIN(Y-E) & Y2=1.1*MAX(Y+E)
XX=FINDGEN(10001)/10000.D0 & XX=0.D0+X1+(X2-X1)*XX
G_ORIG=BETA(0)+BETA(1)*COS(XX)+BETA(2)*SIN(XX)
MITTAA,XX,G_ORIG,AMP_ORIG,MIN_ORIG
RESULT1=AMP_ORIG & RESULT2=MIN_ORIG
; ______________________________________________________________;
FOR QQ=0,ROUNDS-1 DO BEGIN & ERASE ; Bootstrap begins
; ______________________________________________________________; 
!P.TITLE='!6y!Di!N!3'              ; Data and model
SET_XY,X1,X2,Y1,Y2    & SET_VIEWPORT,.1,0.3,0.7,0.9
PLOTERR,T,Y,E,PSYM=4  & PLOT,XX,G_ORIG,PSYM=0
OPLOT,[X1,X2],MAX(G_ORIG)*[1,1],PSYM=0,LINESTYLE=1
OPLOT,[X1,X2],MIN(G_ORIG)*[1,1],PSYM=0,LINESTYLE=1
OPLOT,MIN_ORIG*[1,1],[Y1,Y2],   PSYM=0,LINESTYLE=2
; __________________________________________________________; 
!P.TITLE='!6g!Di!N!3'              ; Model = g
SET_XY,X1,X2,Y1,Y2    & SET_VIEWPORT,0.1,0.3,0.4,0.6
PLOT,T,G,PSYM=4
; __________________________________________________________; 
!P.TITLE='!7e!6!Di!N!3'            ; t and epsilon
Y3=-1.D0*MAX(ABS(EPSILON)) & Y4=-1.D*Y3
SET_XY,X1,X2,Y3,Y4    & SET_VIEWPORT,.1,.3,0.1,0.3
PLOT,T,EPSILON,PSYM=4
; __________________________________________________________; 
!P.TITLE='!7e!6!Di!N!U*!N!3'       ; t and epsilon^*
SET_XY,X1,X2,Y3,Y4    & SET_VIEWPORT,0.4,0.6,0.7,0.9
K=FIX(N*RANDOMU(S,N)) ; Random integers  within [0,K-1]
EPSILON_STAR=EPSILON(K) & E_STAR=E(K)
W_STAR=E_STAR^(-2.0D0)  & Y_STAR=G+EPSILON_STAR ; w^* and y^*
PLOT,T,EPSILON_STAR,PSYM=4
; __________________________________________________________; 
!P.TITLE='!6y!Di!N!U*!N=g!Di!N+!7e!6!Di!N!U*!N!3' ; t and y^*
SET_XY,X1,X2,Y1,Y2    & SET_VIEWPORT,.4,.6,0.4,0.6
PLOTERR,T,Y_STAR,E_STAR,PSYM=4 & Q=DBLARR(3)
NY=CURVEFIT(T,Y_STAR,W_STAR,Q,SIGMAA,FUNCTION_NAME='FU')
G_NOW=Q(0)+Q(1)*COS(XX)+Q(2)*SIN(XX)
PLOT,XX,G_NOW,PSYM=0
MITTAA,XX,G_NOW,AMP_NOW,MIN_NOW
PRINT,'A = ',AMP_NOW,', ... t_min,1 = ',MIN_NOW
RESULT1=[RESULT1,AMP_NOW] & RESULT2=[RESULT2,MIN_NOW]
OPLOT,[X1,X2],MAX(G_NOW)*[1,1],PSYM=0,LINESTYLE=1
OPLOT,[X1,X2],MIN(G_NOW)*[1,1],PSYM=0,LINESTYLE=1
OPLOT,MIN_NOW*[1,1],[Y1,Y2],   PSYM=0,LINESTYLE=2
; ____________________________________________________________; 
!P.TITLE='!6A estimates !3'       ;  A estimates
SET_VIEWPORT,0.7,0.9,0.7,0.9 
Y5=MEAN(RESULT1)-3.*STDEV(RESULT1)
Y6=MEAN(RESULT1)+3.*STDEV(RESULT1)
SET_XY,0.,ROUNDS,Y5,Y6 & IF (QQ GT 1) THEN PLOT,RESULT1,PSYM=4
; _____________________________________________________________; 
!P.TITLE='!6t!Dmin,1!N estimates!3' ; t_min,1 estimates
SET_VIEWPORT,0.7,0.9,0.4,0.6
Y5=MEAN(RESULT2)-3.*STDEV(RESULT2)
Y6=MEAN(RESULT2)+3.*STDEV(RESULT2)
SET_XY,0.,ROUNDS,Y5,Y6 & IF (QQ GT 1) THEN PLOT,RESULT2,PSYM=4
; _____________________________________________________________; 
WAIT,AIKAA                        ; Waits AIKAA seconds
ENDFOR                            ; Bootstrap ends
; _____________________________________________________________; 
PRINT,'     A = ',RESULT1(0),'   +/-',STDEV(RESULT1(1:ROUNDS-1))
PRINT,' t_min = ',RESULT2(0),'   +/-',STDEV(RESULT2(1:ROUNDS-1))
END
; _____________________________________________________________;