!****************************************************************
!***Routine for Cumulative Probabilities for Each Distribution***
!****************************************************************
	SUBROUTINE DIST_CUMPROB(XCALL,PARMDIST,CUMPROB)
	USE SHAREDVARS
	REAL*8 CUMPROB,XCALL,PARMDIST(3)
	REAL*8 XP,X50,XMODE,X0,X100,Z
	REAL*8 SD,SCALE,BETA,LAMBDA,FRACM,KAPPA,MU,SIG,XI,SLP
	REAL*8 PI/3.141592653589793D0/
	REAL*8 NORMZP,NORMPZ
	EXTERNAL NORMZP,NORMPZ
!***
!***Logistic Distribution
!***   Calling Parameter 1 is ECP (XP), P=PTARGET(CLOPT)
!***   Calling Parameter 2 is 1/Scale !!Standard Deviation (SD)
!***   Native Parameter 1 is EC50 (X50)
!***   Native Parameter 2 is Scale Parameter (SCALE) 
!***
	IF(MODELOPT.EQ.1) THEN
          	XP=PARMDIST(1)
		SCALE=PARMDIST(2)
          	X50=XP+SCALE*LOG((1-PTARGET(CLOPT))/PTARGET(CLOPT))
          	CUMPROB=1/(1D0+EXP(-(XCALL-X50)/SCALE))
!***
!***Normal/Gaussian Distribution
!***   Calling Parameter 1 is ECP (XP), P=PTARGET(CLOPT)
!***   Calling Parameter 2 is Standard Deviation (SD)
!***   Native Parameter 1 is EC50 (X50)
!***   Native Parameter 2 is SD 
!***
	ELSEIF(MODELOPT.EQ.2) THEN
		XP=PARMDIST(1)
		SD=PARMDIST(2)
		X50=XP-SD*NORMPZ(PTARGET(CLOPT))
		Z=(XCALL-X50)/SD
		CUMPROB=NORMZP(Z)
!***
!***Rectangular/Uniform Distribution
!***   Calling Parameter 1 is ECP (XP), P=PTARGET(CLOPT)
!***   Calling Parameter 2 is X100-X0 (SCALE)
!***   Native Parameter 1 is X0
!***   Native Parameter 2 is X100 
!***
	ELSEIF(MODELOPT.EQ.3) THEN
		XP=PARMDIST(1)
		SCALE=PARMDIST(2)
		X0=XP-PTARGET(CLOPT)*SCALE
		X100=XP+(1-PTARGET(CLOPT))*SCALE
		IF(XCALL.LE.X0) THEN
			CUMPROB=0D0
		ELSEIF(XCALL.LE.X100) THEN
			CUMPROB=(XCALL-X0)/(X100-X0)
		ELSE
			CUMPROB=1D0
		ENDIF
!***
!***Symmetric Triangular Distribution
!***   Calling Parameter 1 is ECP (XP), P=PTARGET(CLOPT)
!***   Calling Parameter 2 is X100-X0 (SCALE)
!***   Native Parameter 1 is X0
!***   Native Parameter 2 is X100 
!***
	ELSEIF(MODELOPT.EQ.4) THEN
		XP=PARMDIST(1)
		SCALE=PARMDIST(2)
		IF(PTARGET(CLOPT).LE.0.5D0) THEN
			X0=XP-SCALE*SQRT(PTARGET(CLOPT)/2)
			X100=XP+SCALE*(1-SQRT(PTARGET(CLOPT)/2))
		ELSE
			X0=XP-SCALE*(1-SQRT((1-PTARGET(CLOPT))/2))
			X100=XP+SCALE*SQRT((1-PTARGET(CLOPT))/2)
		ENDIF
		IF(XCALL.LE.X0) THEN
			CUMPROB=0D0
		ELSEIF(XCALL.LE.(X0+X100)/2) THEN
			CUMPROB=2*(XCALL-X0)**2/(X100-X0)**2
		ELSEIF(XCALL.LE.X100) THEN
			CUMPROB=1D0-2*(X100-XCALL)**2/(X100-X0)**2
		ELSE
			CUMPROB=1D0
		ENDIF
!***
!***GumbelR (Right Skewed) Distribution
!***   Calling Parameter 1 is ECP (XP), P=PTARGET(CLOPT)
!***   Calling Parameter 2 is scale parameter (BETA) 
!***   Native Parameter 1 is XMODE
!***   Native Parameter 2 is BETA 
!***
	ELSEIF(MODELOPT.EQ.5) THEN
		XP=PARMDIST(1)
		BETA=PARMDIST(2)
		XMODE=XP+BETA*LOG(-LOG(PTARGET(CLOPT)))
		CUMPROB=EXP(-EXP(-(XCALL-XMODE)/BETA))
!***
!***GumbelL (Left Skewed) Distribution
!***   Calling Parameter 1 is ECP (XP), P=PTARGET(CLOPT)
!***   Calling Parameter 2 is scale parameter (BETA)
!***   Native Parameter 1 is XMODE
!***   Native Parameter 2 is BETA 
!***
	ELSEIF(MODELOPT.EQ.6) THEN
		XP=PARMDIST(1)
		BETA=PARMDIST(2)
		XMODE=XP-BETA*LOG(-LOG(1D0-PTARGET(CLOPT)))
		CUMPROB=1-EXP(-EXP((XCALL-XMODE)/BETA))
!***
!***Exponential Distribution with Threshold
!***   Calling Parameter 1 is ECP (XP), P=PTARGET(CLOPT)
!***   Calling Parameter 2 is scale parameter (LAMBDA)
!***   Native Parameter 1 is X0
!***   Native Parameter 2 is LAMBDA 
!***
	ELSEIF(MODELOPT.EQ.7) THEN
		XP=PARMDIST(1)
		LAMBDA=PARMDIST(2)
		X0=XP+LOG(1-PTARGET(CLOPT))/LAMBDA
		IF(XCALL.LE.X0) THEN
			CUMPROB=0D0
		ELSE
			CUMPROB=1-EXP(-LAMBDA*(XCALL-X0))
		ENDIF
!***
!***Asymmetric Triangular Distribution
!***   Calling Parameter 1 is ECP (XP), P=PTARGET(CLOPT)
!***   Calling Parameter 2 is XMODE-X0 (SCALELO)
!***   Calling Parameter 3 is X100-XMODE (SCALEHI
!***   Native Parameter 1 is X0
!***   Native Parameter 2 is X100 
!***   Native Parameter 3 is XMODE
!***
	ELSEIF(MODELOPT.EQ.8) THEN
		XP=PARMDIST(1)
		SCALE=PARMDIST(2)
		FRACM=PARMDIST(3)
		IF(PTARGET(CLOPT).LE.FRACM) THEN
			X0=XP-SQRT(PTARGET(CLOPT)*SCALE**2*FRACM)
			X100=X0+SCALE
			XMODE=X0+FRACM*SCALE
		ELSE
			X100=XP+SQRT((1-PTARGET(CLOPT))*SCALE**2*(1-FRACM))
			X0=X100-SCALE
			XMODE=X0+FRACM*SCALE
		ENDIF
		IF(XCALL.LE.X0) THEN
			CUMPROB=0D0
		ELSEIF(XCALL.LE.XMODE) THEN
			CUMPROB=(XCALL-X0)**2/(X100-X0)/(XMODE-X0)
		ELSEIF(XCALL.LE.X100) THEN
			CUMPROB=1D0-(X100-XCALL)**2/(X100-X0)/(X100-XMODE)
		ELSE
			CUMPROB=1D0
		ENDIF
!***
!***Skewed Logistic Distribution Sastry 
!***   Calling Parameter 1 is ECP (XP), P=PTARGET(CLOPT)
!***   Calling Parameter 2 is Scale Parameter (BETA)
!***   Calling Parameter 3 is Skew Parameter (KAPPA)
!***   Native Parameter 1 is XMODE
!***   Native Parameter 2 is BETA 
!***   Native Parameter 3 is KAPPA
!***
	ELSEIF(MODELOPT.EQ.9) THEN
		XP=PARMDIST(1)
		BETA=PARMDIST(2)
		KAPPA=10D0**PARMDIST(3)
		IF(PTARGET(CLOPT).LE.KAPPA/(KAPPA+1/KAPPA)) THEN
			XMODE=XP+BETA*KAPPA*LOG(2*KAPPA**2/(PTARGET(CLOPT)*(1+KAPPA**2))-1)
		ELSE
			XMODE=XP-BETA/KAPPA*LOG(2*KAPPA**2/(PTARGET(CLOPT)*(1+KAPPA**2)-KAPPA+1)-1)
		ENDIF
		IF(XCALL.LE.XMODE) THEN
			CUMPROB=2*KAPPA**2/(1+KAPPA**2) / (1+EXP(-(XCALL-XMODE)/KAPPA/BETA))
		ELSE
			CUMPROB=KAPPA**2/(1+KAPPA**2) + 2/(1+KAPPA**2) * (1/(1+EXP(-(XCALL-XMODE)*KAPPA/BETA)) - 0.5)
		ENDIF 
!***
!***Skewed Logistic Distribution Aljarrah - POSITIVE SKEW ONLY FOR NOW - SID>0 
!***   Calling Parameter 1 is Location Parameter (Mu)
!***   Calling Parameter 2 is Scale Parameter (Sig)
!***   Calling Parameter 3 is Skew Parameter (Xi)
!***   Native Parameter 1 is Mu
!***   Native Parameter 2 is Sig
!***   Native Parameter 3 is Xi
!***
	ELSEIF(MODELOPT.EQ.10) THEN
		XP=PARMDIST(1)
		SIG=PARMDIST(2)
		XI=10D0**PARMDIST(3)
		MU=XP-SIG*LOG((1-XI*LOG(1-PTARGET(CLOPT)))**(1/XI)-1)
		CUMPROB=1-EXP((1-(1+EXP((XCALL-MU)/SIG))**XI)/XI)
!***
!***End Model Options
!***
	ENDIF
	RETURN
	END SUBROUTINE DIST_CUMPROB


!***
!***Cumulative normal and inverse cumulative normal distributions are based on formula in Sokal and Rohlf, Statistical Tables, 1969
!***

	FUNCTION NORMZP(Z)
      REAL*8 T,B1,B2,B3,B4,B5,C,PI,NORMZP,P,Z
      DATA B1,B2,B3,B4,B5,C,PI/0.319381530D0,-0.356563782D0,1.781477937D0,-1.821255978D0,1.330274429D0,0.2316419D0,0.3989422804D0/
      T=1.0D0/(1.0D0+C*ABS(Z))
      IF (ABS(Z).LT.10) THEN
      	P=1.0D0-PI*EXP(-Z**2/2.0D0)*T*(B1+T*(B2+T*(B3+T*(B4+T*B5))))
      ELSE
      	P=1.0D0
      ENDIF
      IF(Z.LT.0.0) P=1.0D0-P
      NORMZP=P
      RETURN
      END

	FUNCTION NORMPZ(P)
	REAL*8 T,B1,B2,B3,B4,B5,C,PI,NORMPZ,P,PTARG,PP(3),ZZ(3)
	INTEGER IPT,IT
	DATA B1,B2,B3,B4,B5,C,PI/0.319381530D0,-0.356563782D0,1.781477937D0,-1.821255978D0,1.330274429D0,0.2316419D0,0.3989422804D0/
	IF (P.LT.0.0000003) THEN
		NORMPZ=-5.0
		RETURN
	ELSEIF (P.GT.0.9999997) THEN
		NORMPZ=5.0
		RETURN
      ELSEIF(P.GT.0.5D0) THEN
		PTARG=P
	ELSE
		PTARG=1-P
	ENDIF
	ZZ(1)=0.0D0
	ZZ(2)=1.0D0
	ZZ(3)=5.0D0
	DO IPT=1,3      
		T=1.0D0/(1.0D0+C*ZZ(IPT))
		PP(IPT)=1.0D0-PI*EXP(-ZZ(IPT)**2/2.0D0)*T*(B1+T*(B2+T*(B3+T*(B4+T*B5))))
	ENDDO
	DO IT=1,20
            IF(PP(2).GT.PTARG) THEN
			ZZ(3)=ZZ(2)
			PP(3)=PP(2)
		ELSE
			ZZ(1)=ZZ(2)
			PP(1)=PP(2)
		ENDIF
		ZZ(2)=(ZZ(1)+ZZ(3))/2
		T=1.0D0/(1.0D0+C*ZZ(2))
		PP(2)=1.0D0-PI*EXP(-ZZ(2)**2/2.0D0)*T*(B1+T*(B2+T*(B3+T*(B4+T*B5))))
	ENDDO			
	IF(P.GT.0.5D0) THEN
		NORMPZ=ZZ(2)
	ELSE
		NORMPZ=-ZZ(2)
	ENDIF
   	RETURN
	END


